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INTRODUCTION 


The  purpose  of  this  research  is  lo  investigate  the  effects  of  various  electrode 
configurations  on  the  performance  of  the  DPF.  These  studies  focus  on  the  sensitivity  of 
the  current,  density  and  temperature  of  the  plasma  to  different  elecu-ode  geometries.  A 
modularized  version  of  a  DPF  code[l]  previously  used  to  evaluate  the  system  performance 
of  a  D-^He  fueled  device  wiU  be  utilized  for  parametric  studies  in  this  work. 

The  dense  plasma  focus  is  composed  of  a  coaxial  electrode  set  connected  to  a 
high  voltage,  high  current  switching  circuit  (Figure  1).  The  electrodes  consist  of  a 
cylindrical  cathode  surrounding  a  rod  shaped  anode  which  forms  an  annular  gap  between 
the  oppositely  polarized  surfaces.  One  end  of  the  device  is  blocked  off  by  the  fuel  injection 
apparatus  and  insulating  material  needed  to  isolate  the  cathode  from  the  anode.  This  axial 
type  of  coaxial  geometry  is  referred  to  as  the  Mather-type  electrode  configuration,  while 
another  type  of  DPF  electrode  geometry  is  known  as  the  Fillipov  configuration.  The 
Mather  device  accelerates  a  plasma  sheath  axially  down  the  annular  region  of  the  electrodes 
while  a  Fillipov  device  accelerates  a  plasma  sheath  in  a  radial  direction[2].  This  research 
will  be  limited  to  studying  the  Mather  type  of  electrode  configuration  for  space  propulsion 
applications. 

The  DPF's  coaxial  geometry  is  similar  to  other  advanced  propulsion  devices 
like  the  magnetoplasmadynamic(MPD)  and  arcjet  thrusters.  However,  the  MPD  and  arcjet 
do  not  utilize  fusion  power  to  generate  thrust,  and  they  rely  instead  on  expanding  gas  for 
propulsive  power.  The  DPF's  dependence  on  fusion  power  requires  operational  currents 
that  are  far  in  excess  of  the  current  required  by  the  MPD  or  arcjet.  Another  difference  is 
that  the  DPF  is  not  operated  in  a  steady  state  mode,  rather  it  is  a  pulsed  power  device. 

The  principal  operation  of  the  DPF  is  fairly  straightforward;  fusion  fuel  is  first 
injected  into  the  annulus  between  both  electrodes  and  an  arc  is  established  across  the  fuel 
filled  gap.  This  arc  ionizes  the  fuel  and  forms  a  plasma  sheath  which  propagates  down  the 
annulus  and  runs  out  the  end  of  the  device.  As  the  outer  rim  of  the  sheath  detaches  from 
the  cylindrical  cathode,  magnetohydrodynamic  instabilities  cause  the  sheath  to  collapse  in 
on  itself  thereby  forming  a  small,  hot,  highly  dense  volume  of  plasma.  It  is  in  this  region 
that  the  fusion  reactions  occur,  which  liberate  the  necessary  energy  used  for  thrust. 
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Figure  I 

Diagram  of  Dense  Plasma  Focus  Device 
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The  operational  cycle  of  the  DPF  can  be  broken  down  into  several  distinct 
phases  as  illustrated  in  Figure  2. 

1.  Breakdown  phase  -  Gaseous  fuel  is  injected  into  the  annular  region  prior  to  arc  initiation. 
Capacitor  bank  is  discharged  across  electrodes  and  initiates  symmetrical  arc  between 
cathode  cylinder  and  anode  bar.  Fill  gas  is  ionized  during  breakdown  and  plasma  sheath 
is  formed. 

2.  Rundown  phase  -  Arc  current  induces  azimuthal  magnetic  field  Be  around  the  anode  bar. 
The  J  X  Be  force  accelerates  plasma  sheath  down  the  length  of  the  anode.  A  fraction  of 
the  fill  gas  is  entrained  by  the  propagating  plasma  sheath. 

3.  Pinch  phase  -  Plasma  sheath  reaches  the  end  of  the  electrodes  and  the  J  x  Be  force 
initiates  the  radial  compression  of  the  sheath.  Collapsing  sheath  focuses  towards  the 
central  axis  of  the  anode  forming  a  hot,  high  density  plasma  where  fusion  reactions  are 
to  take  place. 

After  the  pinch  is  formed,  it  is  vulnerable  to  various  types  of  plasma  instabilities 
which  will  distort  and  eventually  disrupt  the  pinch  region  due  to  MHD  instabilities.  These 
instabilities  include  both  the  m  =  0  "sausage"  and  m  =  I  "kink"  instabilities.  The 
disruption  of  the  pinch  region  marks  the  end  of  one  cycle  of  the  DPF  operation,  and  this 
cycle  is  designed  to  be  repeated  providing  a  "continuous"  mode  of  firing. 

Advanced  propulsion  technology  development  began  in  the  I950's  with  the 
development  of  nuclear  and  ion  propulsion  systems.  These  concepts  came  under  scrutiny 
because  of  the  potential  advantages  over  chemical  rocket  concepts[3,4].  Advanced 
concepts  relying  on  nuclear,  ion,  and  magnetoplasmadynamic  schemes  could  provide  the 
increase  in  thrust  and  specific  impulse  necessary  for  manned  interplanetary  travel. 

Nuclear  propulsion  offers  much  more  energy  per  unit  mass  than  conventional 
chemical  rocket  concepts.  Chemical  rockets  are  limited  by  the  chemical  bond  energy  of  the 
fuel,  whereas  nuclear  rockets  rely  on  fission  energy.  This  available  nuclear  energy  far 
exceeds  that  produced  by  even  the  most  energetic  chemical  rocket  fuels.  Project  Rover[5] 
was  one  of  the  first  studies  to  examine  the  feasibility  of  nuclear  rocket  propulsion.  Thrust 
for  the  nuclear  rocket  is  produced  by  heating  propellant  in  a  heat-exchanger  reactor  and 
exhausting  the  hot  gas  through  a  nozzle.  Generally,  hydrogen  is  the  propellant  of  choice 
because  of  it's  low  molecular  weight.  The  rocket  itself  consists  of  a  solid  core  reactor  with 
numerous  flow  channels  to  accommodate  the  passage  of  propellant  thicugh  the  core.  The 
Rover  project  was  concerned  with  producing  roughly  twice  the  specific  impulse  of  the  most 
efficient  chemical  rocket,  as  well  as  providing  greater  thrust  for  an  increased  payload. 
Open  air  tests  were  conducted  with  the  Kiwi-A  device  in  1959-60,  and  it  concluded  with 
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the  recommendation  of  further  study  into  nuciear  propulsion.  Another  nuclear  rocket 
project  that  followed  in  the  footsteps  of  the  Rover  project  was  the  NERVA  project.  The 
goal  of  the  NERVA  project  was  to  produce  a  solid  core  powered  nuclear  propulsion  system 
that  would  be  capable  of  the  high  specific  impulse  needed  for  interplanetary  travel.  The 
NERVA  research  began  in  the  i960's,  but  fell  victim  to  funding  cuts  in  the  early  1970's. 

Ion  propulsion  is  another  advanced  concept  that  was  theoretically  predicted  to 
be  superior  to  chemical  propulsion.  Ion  rockets  rely  on  the  electrostatic  acceleration  of  ions 
for  thrust[4].  An  ion  rocket  is  a  form  of  electric  propulsion  and  therefo  requires  a  power 
source  to  produce  ions  and  generate  the  electrostatic  field  used  to  accelerate  them.  Ion 
propulsion  performance  exceeds  that  of  the  nuclear  rocket  in  terms  of  travel  time  and 
specific  impulse,  but  it  lacks  the  thrust  to  weight  capacity  of  nuclear  and  chemical  systems. 
Consequently,  an  ion  propulsion  system  must  first  be  placed  in  orbit  before  it  can  engaged. 

Magneto-plasmadynamic(MPD)  thrusters  are  a  form  of  electrical  propulsion  that 
can  provide  performance  similar  to  ion  powered  devices.  MPD  type  thrusters  utilize  a 
current  and  magnetic  field  to  produce  an  electromagnetic(JxB)  force  that  is  used  to 
accelerate  a  highly  ionized  fuel.  These  thrusters  basically  consist  of  an  cathode  and  anode, 
either  in  a  parallel  rail  or  coaxial  electrode  configuration,  and  a  power  supply.  Operation  of 
this  type  of  thruster  is  similar  to  that  of  the  DPF  which  is  described  in  the  previous  section. 
One  of  the  major  differences  between  the  MPD  and  the  DPF  is  that  the  DPF  electrode 
polarity  is  reversed  from  the  MPD  electrode  polarity.  The  arcjet  which  is  currently  under 
research  at  the  NASA  Lewis  Center  belongs  to  this  family  of  electric  thruster. 

One  alternative  concept  that  has  not  been  examined  as  extensively  as  the 
previously  mentioned  schemes  is  fusion  propulsion.  Fusion  propulsion  relies  on  the 
enormous  amounts  of  energy  produced  in  thermonuclear  reactions  to  produce  thnist. 
Theoretical  studies  have  been  conducted  on  different  fusion  propulsion  schemes  based  on 
inertial  and  magnetic  confinement  concepts[6,7].  The  Air  Force  has  undertaken  a  study 
involving  the  dense  plasma  focus  device  and  its  application  to  space  propulsion.  The 
feasibility  of  designing  and  operating  such  a  system  is  being  researched  at  the  Phillips 
Laboratory  in  New  Mexico  and  Purdui  University.  The  adaptation  of  the  den.se  plasma 
focus  device  for  space  propulsion  involves  an  increased  current  and  material  requirement 
that  is  essential  for  operation.  Operating  currents  on  the  order  of  tens  of  mega  amperes 
must  be  applied  in  order  to  achieve  a  fusion  reaction  with  a  D-^He  fueled  device.  The 
current  capability  of  the  source  circuit  as  well  as  the  durability  of  electrode  and  insulator 
material  represent  major  engineering  limitations.  The  physics  oi  the  pinch  dynamics,  and 
the  scaling  laws  relating  plasma  temperature  and  capacitor  ma.ss  to  external  current  are 
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additional  problems  which  require  further  research,  but  are  beyond  the  scope  of  the  curreni 
study. 

Several  plasma  focu.s  devices  have  been  constructed  and  tested  at  various 
research  institutions  around  the  worldlH).  The  Livermore- 1  e.xpenment  which  forms  the 
basis  of  this  study  is  a  1.2  MA  device  that  was  used  to  conduct .  ludies  into  the  rundown, 
collapse  and  pinch  phases  of  operation.  Other  devices  include  the  Frascati  plasma  fcx'us  m 
Italy  which  operated  at  a  peak  curreni  of  2.8  MA,  and  the  Poseidon  plasma  focus  device  in 
Stuttgart  Germany  which  reached  a  maximum  curreni  of  4.9  MA19].  These  devices  u.sed 
annealed  copper  electrodes  and  pyrex  or  ceramic  insulators.  These  maieriaLs  help  determine 
a  failure  limit  for  the  current  generation  of  devices.  In  order  to  achieve  ignition,  a  D-^He 
fueled  device  should  be  supplied  with  a  total  mput  current  in  the  10-20  MA  range[  1].  Tins 
limit  far  e.xceeds  the  capacity  of  any  exLsting  pla.sma  focus  device.  However,  the  SHIVA 
implosion  experiment  at  the  Air  Force  Weapons  research  lab[  10]  reached  a  peak  operating 
current  of  approximately  12.. 3  MA.  This  experimeni  establishes  a  precedent  for  very  high 
current  capability  for  a  pulsed  power  system  application.  The  SHIVA  experiment  relied  on 
a  120  kV,  9.4-MJ  capacitor  bank  to  supply  the  necessary  curreni  to  drive  a  dynamic  coaxial 
vacuum  inductive  store.  This  system  was  u.sed  to  supply  a  current  pulse  to  a  cylindrical 
implosion  load. 

The  dense  plasma  focus  has  been  the  subject  of  numerous  studies  since  its 
inception  in  the  1960's.  and  the  majority  of  the.se  studies  have  been  conducted  on  the 
plasma  pinch  phenomena  and  neutron  production  mechanisms!  11,12).  Parametric  studies 
have  been  conducted  with  the  goal  of  optimizing  the  circuit  parameters  and  filling  pressure 
to  maximize  the  temperature  and  density  of  the  pinch  region!  13].  Variation  of  electrode 
configuration  and  the  consequent  effect  on  the  pinch  have  not  been  as  extensively  studied 
as  other  plasma  focus  phenomena. 

The  electrode  configuration  plays  a  key  role  in  the  formation  of  the  plasma 
pinch  region.  Geometry  and  .sizing  of  the  electrodes  affect  the  density,  dynamic 
inductance,  and  therefore  the  current  of  the  propagaung  plasma  sheath.  The  characteristics 
of  the  current  pulse  can  be  used  to  determine  the  electrode  length  in  order  to  ensure  that 
maximum  current  is  reached  at  the  end  of  the  acceleration  phase!  13J.  Expenmentatior  with 
various  plasma  focus  devices  such  as  the  Livermore-I  experiment  has  yielded  information 
about  the  current  history  of  the  propagating  plasma  sheath.  Subsequent  studies  have  dealt 
with  the  presence  of  a  leakage  current  that  occurs  over  the  surface  of  the  insulator  during 
operation[8].  This  leakage  current  degrades  the  performance  of  the  focus  for  operating 
currents  in  the  MA  range.  The  high  operating  current  range  al,so  has  an  effect  on  the 
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insulator  surface  over  which  the  initial  breakdown  arc  occurs(14].  In  a  22  kV,  1-kJ 
Mather-type  plasma  focus,  the  pyrex  insulator  used  suffered  surface  alterations  and  metallic 
deposition  due  to  the  temperature  in  the  plasma  sheath,  llie  device  utilized  brass  electrodes 
and  reached  a  peak  current  of  approximately  0.1  MA.  After  successive  firings  of  the 
device,  varying  stages  of  insulator  erosion  were  observed.  This  degradation  is  inuicative  of 
the  susceptibility  of  the  insulator  to  very  high  operating  currents  and  is  a  key  factor  in 
determining  the  limiting  performance  of  upscaled  models  of  very  high  current  devices. 

Various  models  have  been  used  to  predict  current  histories  and  sheath  velocities 
in  the  rundown  phase  of  operation.  Different  studies  have  utilized  various  modeling 
techniques  for  the  plasma  focus  rundown,  including  an  equivalent  circuit 
representation[8,15],  and  two-dimensional  MHD  calculations.  The  equivalent  circuit 
models  provide  a  simpler  model  which  couples  electrical  circuit  equations  to  dynamic 
sheath  parameters  while  the  MHD  codes  use  a  two-Ouid  model  with  two  dimensional 
effects  included[16,17].  Both  methods  can  be  used  to  predict  the  behavior  of  the  rundown 
phase  up  to  the  collapse  region,  but  MHD  theory  breaks  down  after  the  collapse  since  it 
cannot  predict  the  pinch  formation,  dimension  and  the  current  behavior  in  the  pinch  region. 

This  work  is  being  conducted  in  order  to  establish  the  feasibiliiy  of  the  dense 
plasma  focus  as  a  viable  space  propulsion  concept.  The  enormous  potential  of  fusion 
power  in  this  application  is  explored  for  the  dense  plasma  focus  device.  To  achieve  a  good 
fusion  bum  during  the  pinch  phase,  currents  on  the  order  of  20  MA  must  be  supplied  to  the 
device.  Present  plasma  focus  experiments  deal  with  sheath  currents  that  are  only  in  the  I  - 
4  MA  range.  Therefore,  a  realistic  model  must  be  developed  for  the  upscaled  parameters  of 
a  dense  plasma  focus  system  that  operates  in  the  very  high  current  regime. 

In  order  to  accurately  account  for  the  physics  of  the  breakdown  and  the 
mndown  phases  of  operation,  a  or.e-dimensional  transient  simulator  code  will  be  developed 
for  integration  with  a  previously  developed  code  for  the  pinch  phase  of  the  plasma 
focus[I].  This  model  will  be  based  on  an  equivalent  circuit  representation  of  the  dense 
plasma  focus  device.  Various  features  will  be  incorporated  into  the  code  to  account  for 
certain  phenomena  that  have  been  experimentally  observed.  A  leakage  current  branch  in  the 
equivalent  circuit  is  used  to  provide  a  realistic  loss  component  to  the  transient  code.  The 
current  of  the  plasma  sheath  prior  to  the  pinch  phase  ideally  should  be  at  a  maximum  in 
order  to  optimize  conditions  for  a  fusion  reaction  to  take  place.  In  order  to  reduce  the 
current  damping  during  operation,  a  crowbar  switch  was  also  added  to  the  equivalent 
circuit.  The  equivalent  circuit  is  reduced  to  a  system  of  equations  which  is  solved  using 
LU  decomposidon.  An  objective  will  be  to  link  the  electrical  performance  of  the  device  to 
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the  plasma  modeling  of  the  sheath.  The  Snowplow  model  is  utilized  to  calculate  rundown 
velocity  and  the  mass  of  fill  gas  entrained  in  the  propagating  sheath.  The  dynamic  sheath 
inductance  and  resistance  predicted  by  the  Snowplow  model  is  coupled  to  the  equivalent 
circuit  system  of  equations.  This  approach  provides  a  realistic  means  of  predicting  plasma 
focus  performance  in  the  20  MA  current  range.  Calculations  will  then  be  carried  out  on 
different  electrode  geometries  with  the  goal  of  obtaining  a  more  accurate  representation  of 
the  performance  envelope  of  the  dense  plasma  focus. 


MODELING  OF  THE  DENSE  PLASMA  FOCUS  (DPF) 


Equivalent  Circuit 

The  dense  plasma  focus  device  can  be  modeled  using  an  equivalent  circuit 
representation.  The  model  shown  in  Figure  3  displays  the  circuit  parameters  for  the 
discharge  circuit  and  the  plasma  sheath.  The  values  C,  Lq  and  Rq  are  the  external 
capacitance,  inductance  and  resistance,  respectively.  The  dynamic  inductance  and 
resistance  of  the  plasma  sheath  are  represented  by  Lpp  and  Rpp. 


External  Discharge  Plasma  Sheath 

Circuit  Parameters  Circuit  Parameters 


Figure  3 

Equivalent  Circuit  Model  for  Dense  Plasma  Focus 


The  resistance  Rl  that  is  included  in  the  sheath  parameters  depicts  the  leakage  current 
around  the  insulator  that  occurs  during  the  rundown  phase.  The  leakage  current  results 
from  the  formation  of  an  arc  across  the  insulator  surface,  and  is  responsible  for  reducing 
the  actual  current  delivered  to  the  plasma  sheath.  Rcr  and  Lcr  are  the  components  of  the 
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crowbar  switch  of  the  circuit.  The  crowbar  switch  is  designed  to  reduce  the  fust  damping 
of  the  current  that  occurs  during  the  capacitor  discharge. 

At  time  t  =  0  ■ ,  switches  S I  and  S2  are  open,  then  S 1  is  closed  at  t  =  0  and  the 
capacitor  begins  discharging.  When  the  current  reaches  its  maximum  value  at  t  =  S2 
is  closed  to  reduce  the  damping  effect.  The  values  of  the  circuit  parameters  can  be  chosen 
to  provide  either  an  underdamped,  overdamped,  or  critically  damped  transient  response. 
For  this  work,  the  external  circuit  parameters  will  remain  fixed  while  the  crowbar  values 
and  the  sheath  values  will  be  optimized  in  order  to  provide  maximum  current  delivery  to  the 
sheath. 


Companion  Circuit  Model 

The  modeling  of  the  reactive  elements  in  the  circuit  is  done  using  the  companion 
circuit  model.  In  the  companion  circuit  model,  both  the  inductors  and  capacitors  in  a  circuit 
can  be  reduced  to  an  equivalent  or  "companion"  representation.  This  is  accomplished  by 
making  use  of  the  basic  voltage  and  current  relations  for  the  inductor  and  capacitor. 


V  =  L^ 
dt 


(1) 


and 


l  =  C^ 

dt 


(2) 


Both  relations  are  evaluated  as  integrals  with  respect  to  time. 


The  solution  to  these  integral  equations  can  be  obtained  by  using  a  trapezoidal 
approximation  for  either  the  voltage  or  current  waveform  as  illustrated  in  Figure  4. 
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Figure  4 

Trapezoidal  Approximation  of  Waveform 


The  value  of  the  function  F(n)  can  be  ascertained  by  summing  the  values  of  each 
trapezoidal  panel  over  a  specified  number  of  time  steps(n).  Equation  5  denotes  the  value  of 
a  single  trapezoidal  panel  for  a  single  time  step  h. 

AI  =  0.5h{F(n)  +  F(n+l}).  (5) 

Usage  of  the  trapezoidal  rule  allows  the  inductor  and  capacitor  to  be  reduced  to  the 
companion  models[18]  shown  in  Figure  5. 


ReacCive 

Component 


Companion 

Equivalent 


Figure  5 

Equivalent  Companion  Circuits 
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These  models  form  the  basis  of  the  governing  equations  for  the  inductor  and  capacitor, 
with  the  (n+1)  terms  representing  the  new  time  values  and  the  (n)  terms  representing  the 
old  time  forcing  function  values. 


.  I(n+1)  =  V(n+1)  +  ^  V(n)  +  I(n)  . 

and  2L  2L 


V(n+l)  =  ^I(n+l)+  ^I(n)  +  V{n) 


(6) 

(7) 


Ctevelopment  of  Transient  Circuit  Code 

Development  of  a  transient  code  was  deemed  necessary  in  order  to  model  the 
electrical  and  plasma  dynamic  characteristics  of  the  dense  plasma  focus.  In  the  transient 
code,  the  dynamic  physics  of  the  sheath  are  coupled  with  the  plasma  inductance  and  current 
response  of  the  equivalent  circuit.  This  coupling  allows  the  device  to  be  modeled  by  a 
modified  transient  circuit  solver. 

The  equivalent  circuit  can  be  reduced  to  a  system  of  equations  represented  by 

1  =  YV.  (8) 

where  I  is  the  current  vector,  Y  is  the  conductance  matrix  and  V  is  the  voltage  vector.  The 
Y  matrix  is  formed  using  the  conductances  of  the  resistors  and  the  trapezoidal 
approximation  to  model  the  reactive  elements  of  the  circuit.  Once  the  Y  matrix  is  formed, 
the  initial  conditions  are  input  into  the  I  vector,  then  the  system  is  solved  using  LU 
decomposition.  The  Y  matrix  and  I  vector  are  updated  every  time  step  to  account  for  the 
dynamically  changing  inductance  and  resistance  of  the  sheath.  Each  new  time  current  is 
used  to  calculate  a  new  rundown  velocity  as  well  as  a  new  sheath  inductance,  resistance 
and  temperature.  This  calculation  is  marched  through  time  until  the  sheath  reaches  the  end 
of  the  device,  then  the  end  values  of  the  current,  rundown  velocity  and  plasma  temperature 
are  passed  to  another  code  which  calculates  the  pinch  phase  dynamics. 

The  operational  cycle  of  the  rundown  phase  can  be  described  by  two  distinct 
circuits,  one  without  the  crowbar  for  0  <  t  <  tn,ax(Figure  6),  and  the  other  with  the  crowbar 
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switched  in  for  lime  t  >  tma)i(Pigure  7)  These  two  circuits  will  yield  conductance  mauices 
of  differing  sizes  and  it  is  desirable  to  keep  them  segregated  to  allow  for  ease  of  calculation. 


Figure  6 

Equivalent  Circuit  for  0  <  t  <  t^ax 


The  circuit  for  0  <  i  <  imax  is  reduced  to  a  system  of  equations  by  summing  the 
currents  into  each  node. 


(9) 

Ko 

(10) 

Ic  +  Ilo  =  0  , 

(11) 

and 

-fc+j^  =  0. 

KpF 

(12) 

This  system  is  augmented  by  an  additional  number  of  equations,  each  representing  the 
companion  circuit  model  of  a  reactive  element 

Idn)  +  ^  V3(n)  =  ^  V3{n+ 1 )  -  Ic{n+ 1 ) . 
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-  rf  -  (V3{n)  -V.(n))  -  Mn)  =  { V3(n+ 1 )  -  V,{n+ 1 ))  -  Mn+ 1 ) . 

and  2Lo  2Lo  (14) 

The  resultant  matrix  that  is  formed  by  both  these  sets  of  equations  is  used  for  calculating 
the  new  time  currents  and  voltages  in  the  circuit  by  LU  decomposition. 

Once  the  crowbar  branch  of  the  circuit  is  switched  in  at  t  >  t^ax-  it  new  conductance 
matrix  must  be  formed  to  account  for  the  modified  circuit.  The  equivalent  circuit  with  the 
crowbar  switched  in  is  displayed  in  Figure  7. 


Lpf 

V4 

Rpf 


Figure  7 

Equivalent  Circuit  for  i  >  tmax 


Switching  the  crowbar  into  the  circuit  introduces  one  extra  node  and  2  extra  elements  to  the 
basic  equivalent  circuit.  Two  extra  nodal  equations  and  the  additional  crowbar  current  are 
added  to  the  previous  set  formulated  for  the  basic  equivalent  circuit  without  the  crowbar. 


•iu.  +  ^  =  o. 


(15) 


1l» (“+‘>  =  3?- ( V,  (n+l )  -  Vs  (n+ 1 ))  +  ^ ( V,  (n)  -  V,  (n))  +  (n) . 

and  (16) 

These  additional  nodal  equations  allow  the  previous  7x7  syvStem  (Figure  8)  of  equations  to 
be  increased  to  a  9x9  system  (Figure  9). 
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These  matrices  were  integrated  into  a  transient  code  developed  specifically  for  the 
rundown  phase  of  the  dense  plasma  focus.  Initial  conditions  are  read  from  an  input  deck 
and  the  appropriate  arrays  are  initialized.  The  code  is  first  run  through  the  problem  with  the 
mode  flag  set  equal  to  3.  This  indicates  that  the  code  is  checking  for  the  time  that  peak 
pinch  current  is  reached  in  the  equivalent  circuit  without  the  crowbar.  Once  the  code  cycles 
through  this  check  mode,  the  time  that  peak  current  is  reached  {tnm)  is  stored  and  the  flag 
is  reset  to  1  and  the  problem  is  restarted.  The  flag  remains  set  at  1  until  tmax  is  reached, 
then  the  crowbar  is  switched  in  and  the  flag  is  set  to  2.  The  calculation  is  now  marched 
through  time  until  the  propagating  plasma  sheath  reaches  the  end  of  the  anode. 

This  code  tracks  the  electrical  and  plasma  parameters  of  the  device  through  the 
duration  of  the  rundown  phase.  As  mentioned  earlier,  the  7x7  matrix  is  used  before  the 
plasma  current  has  reached  its  peak,  and  the  9x9  matrix  is  used  after  peak  current  has  been 
reached  and  the  crowbar  has  been  switched  in.  During  each  time  step  the  matrix  is 
initialized  and  solved  by  LU  decomposition,  and  the  new  time  values  are  used  to  calculate  a 
new  sheath  inductance,  rundown  velocity,  plasma  temperature  and  plasma  resisumce. 
These  values  are  recycled  back  to  initialize  the  matrix  for  the  next  time  step  and  the  process 
repeats  itself  until  the  sheath  runs  off  the  tip  of  the  anode.  A  flow  chart  depicting  the  order 
of  operation  in  the  transient  code  is  displayed  in  Figure  10. 

Initial  Testing  of  Transient  Circuit  Code 

The  validity  of  the  transient  LU  solver  was  established  by  running  a  test  case 
for  a  circuit  with  static  parameters  and  contrasting  the  calculated  results  with  the  analytical 
solution  and  output  from  SPICE.  The  SPICE  (Simulation  Program  with  Integrated  Circuit 
Emphasis)  program[19]  is  a  circuit  simulator  that  was  developed  for  various  types  of  circuit 
analysis,  including  linear  ac,  nonlinear  transient  and  nonlinear  dc  analysis. 

The  test  circuit  was  patterned  after  the  equivalent  circuit  displayed  in  Figure  6. 
Circuit  parameter  values  were  set  to  provide  as  close  a  resemblance  to  the  Livermore  I 
parameters  as  possible.  Table  1  shows  the  vaiious  parameters  for  Livermore  I  and  the 
static  lest  circuit  used. 

Table  1.  Comparison  of  Parameters  between  Static  Testing  Circuit  and  the  Livermore-I 

Experiment 


DEVICE 

C(F) 

Vo(V) 

Lo  (H) 

Rl  («) 

rpf  m 

LpF  (H) 

Livermore-I 

3.55E-4 

27000 

.005 

2.5E-8 

.12 

3.55E-4 

27000 

.005 

2.5E-8 

.12 

l.OE-7 

l.OE-8 
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The  7x7  System  of  Equations  for  the  Equivalent  Circuit  at  0  <  t  <  t, 
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Figure  9 

The  9x9  System  of  Equations  for  the  Equivalent  Circuit  at  t  >  t. 


The  results  from  the  test  runs  on  the  static  circuit  are  displayed  m  the  following 
graphs.  Figure  1 1  shows  the  dual  plot  of  ioad  current  for  both  SPICE  and  LU  calculated 
solutions.  It  is  apparent  that  the  overlay  of  both  plots  make  them  virtually  mdistingui.shable 
from  each  other.  The  relative  error  between  the  two  calculations  is  shown  in  Figure  12. 
There  is  excellent  agreement  between  SPICE  and  the  LU  calculated  solutions  (note  that  the 
percent  error  is  normalized  for  the  SPICE  calculation).  The  oscillating  nature  of  the  relative 
error  is  tliought  to  be  due  to  perturbations  in  the  SPICE  calculated  result.  The  next  graph 
(Figure  13)  shows  the  relative  error  between  the  LU  caLulation  and  an  analytical  solution 
that  was  derived  using  Laplace  Transformations  and  Kirchoffs  Voltage  Law.  Again  there 
is  excellent  agreement  between  the  LU  calculation  and  the  analytical  solution  and  there  is  no 
evidence  of  the  oscillation  present  in  the  Figure  12.  An  actual  graph  between  the  analytical 
value  and  the  LU  solution  was  omitted  since  the  curves  cannot  be  distinguished  from  each 
other. 

The  crowbar  switch  was  a  feature  that  was  added  to  the  transient  code  in  order  to 
reduce  the  d.‘’.nping  of  current  in  the  response.  The  switching  in  of  the  crowbar  branch 
increases  the  effective  decay  constant  of  the  current  response,  thus  keeping  the  current  as 
high  as  possible  for  as  long  as  possible.  The  effectiveness  of  the  crowbar  is  determined  by 
the  choice  of  crowbar  inductance  and  resistance.  The  maximum  effect  can  be  achieved 
when  the  crowbar  parameters  are  of  the  same  order  of  magnitude  as  the  sheath  parameters. 
A  test  case  with  the  crowbar  switched  in  was  run  and  contrasted  with  the  results  from  a 
case  run  without  the  crowbar.  Figure  14  displays  both  plots  and  it  should  be  noted  that  the 
curve  for  the  crowbar  is  for  the  median  case  when  the  crowbar  and  load  impedances  are 
identical. 

The  effect  of  the  crowbar  that  is  evident  in  Figure  14  suggests  that  the  addition  of  a 
crowbar  switch  would  be  beneficial  to  the  maximization  of  current  delivered  to  a  load. 
There  are  however  some  engineering  concerns  that  may  limit  the  usage  of  the  crowbar  for 
plasma  driven  devices.  One  concern  is  that  in  order  to  maximize  the  effect  of  the  crowbar, 
the  crowbar  resistance  and  inductance  must  be  of  the  same  order  of  magnitude  as  that  of  the 
load.  In  devices  that  utilize  a  body  of  plasma  as  the  load,  the  resistance  and  inductance  are 
typically  very  small  values,  and  it  may  be  hard  to  produce  an  adequate  crowbar  switch 
which  would  significantly  decrease  damping.  This  switch  would  also  have  to  be  durable 
enough  to  withstand  flow  currents  in  the  10^  ampere  range.  Assuming  that  these  problems 
can  be  resolved,  the  crowbar  switch  could  prove  to  be  very  beneficial  to  peak  current 
stabilization  in  the  dense  plasma  focus. 
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Figure  1 1 

Cunreni  Plot  for  Test  Circuit  using  SPICE  and  LU  Solver 


Figure  12 

Percent  Error  between  SPICE  and  LU  Solutions 


Figure  13 

Percent  Error  between  LU  and  Analytical  Solutions 


Figure  14 

Calculated  Load  Currents  for  Circuit  with,  and  without  Crowbar. 


21 


Plasma  Dynamics 

The  dynamic  quantities  associated  with  the  propagating  plasma  sheath  in  the  DPF 
device  are  the  sheath  inductance,  resistance,  temperature  and  rundown  velocity.  Relations 
for  these  quantities  allow  die  dynamics  of  the  plasma  sheath  to  be  coupled  to  the  equivalent 
circuit  model.  The  snowplow  model  [20]  is  used  to  calculate  rundown  velocity  in  the 
transient  code.  This  model  is  used  to  account  for  the  mass  entrainment  of  the  fill  gas  as  the 
plasma  sheath  propagates  down  the  anode.  The  simple  MHD  model  does  not  take  this 
entrainment  effect  into  account  and  is  perhaps  too  ideal  for  a  realistic  calculation. 

The  snowplow  model  assumes  that  the  plasma  sheath  is  an  impermeable  surface 
that  absorbs  fill  gas  as  it  propagates  down  the  electrode  annulus.  Any  mass  that  the  arc 
surface  encounters  as  it  sweeps  down  the  anode  is  absorbed  into  the  sheath.  The 
snowplow  model  can  be  expressed  using  Newton’s  law  of  motion  in  which  the  time  rate 
change  of  momentum  is  equal  to  the  sum  of  forces  acting  on  a  body.  The  limitation  of  the 
snowplow  model  is  that  it  assumes  total  mass  entrainment  which  does  not  totally  reflect  the 
physical  nature  of  the  sheath  during  rundown.  One  expects  that  while  the  majority  of  the 
fill  gas  is  swept  into  the  sheath,  there  is  a  certain  fraction  that  is  not  entrained  during 
rundown.  This  difference  does  not  provide  a  significant  discrepancy  in  the  calculated 
results  for  the  Livermore-I  plasma  focus  test  case  that  is  shown  later  in  this  text.  Starting 
with  the  general  momentum  equation: 


3pv 

~ 


+  V-(p  V  v)  =  -Vp  -  V  t  +  pg, 


(17) 


where  p  is  the  initial  fill  gas  density  and  v  is  the  velocity.  The  waU  shear  and  body  force 
terms  are  neglected  under  the  assumption  that  the  magnetic  pressure  is  the  primary  driving 
force.  By  taking  only  the  axial  components  and  integrating  over  the  constant  volume  of  the 
sheath  during  rundown,  a  general  expression  for  the  momentum  balance  can  be  obtained 

=  (18, 


where  the  left  hand  side  of  the  equation  is  the  total  derivative  representation  of  the  time 
dependent  and  convective  momentum  terras.  The  total  force  term  on  the  right  hand  side  of 
Equation  (18)  includes  magnetic  force  on  the  sheath,  particle  pressure  of  the  plasma  and  the 
frictional  force  on  the  plasma  mass.  The  frictional  force  is  negligible  in  comparison  to  the 
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magnetic  force  and  particle  pressure.  The  particle  pressure  will  be  treated  later  in  the  text 
when  the  liftoff  current  is  calculated.  The  magnetic  force  is  expressed  by  the  equation  of 
the  electromagnetic  force  density  integrated  over  the  volume  of  the  sheath, 

Fmag  ~  I  j  X  B  dX  . 

Jvoi  (19) 

The  work  done  by  the  propagating  plasma  sheath  as  it  advances  a  distance  z  down  the 
anode  is  expressed  by  the  next  equation. 


w=  1  F^agdz. 

Jo 

(20) 

Assuming  that  all  of  the  inductive  energy  (including  energy  stored  in 

the  field)  of  the 

plasma  sheath  is  going  into  driving  the  plasma,  the  following  expressions  can  be  obtained: 

Wind=^Ll2. 

(21) 

W  =  Wind, 

(22) 

F  _  1  L  l2 

t^mag  -  2  Z  ■ 

(23) 

Expanding  the  left  hand  side  of  Newton's  law  (Equation  (18))  and  substituting  the  force 
term  with  Equation  (21)  gives  the  following  relation. 

mz  +  rhz= 

2  z  (24) 

In  this  expression  the  first  term  on  the  left  hand  side  of  the  equation  is  the  conventional 
acceleration  term,  while  the  second  term  designates  the  mass  accumulation  of  the  sheath. 

The  plasma  sheath  initially  has  some  small  mass  mo  associated  with  it  at  the 
beginning  of  the  breakdown  phase.  The  initial  mass  quantity  is  dependent  on  the  sheath 
thickness  and  length  at  the  initiation  of  the  arc  sheet.  Once  the  sheath  starts  propagating 
down  the  device,  the  mass  accumulation  rate  is  determined  by  the  rundown  velocity,  gas 
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density  and  the  dimensions  of  the  annular  channel.  An  assumption  is  made  that  the  gas 
density  is  uniform  throughout  the  transient.  The  total  mass  term  m  can  be  expressed  as 
follows; 


m  =  mo  +  Area  I  p  ^  dt 


{ 


dt 


where  Area  —  Jt  {r|  -  rj) 


(25) 

(26) 


The  area  is  calculated  for  the  annular  dimensions  of  the  electrodes  with  the  final  mass 
expression  with  p  being  constant 

m  =  mo  +  Tc(r2- r|)p  z<t)  (27) 

Newton's  law  can  be  also  be  expressed  in  the  following  reduced  form 

•ifmltlzlt)]  =  I^(t} 

dt^^'^'^  2  z  ''  (28) 

Integrating  with  respect  to  time  gives 


Mt)I^(t)dt 

z(t) 


(29) 


Now  we  have  an  expression  for  the  rundown  velocity  which  accounts  for  mass  entrainment 
in  the  sheath. 


z  = 


2  m 


z(t) 


(30) 
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This  expression  is  implemented  numerically  in  the  transient  code  using  the  U'apezoidal  rule 
to  approximate  the  product  of  L  and  P  at  each  discrete  time  step.  This  gives  a  new  time 
rundown  velocity  which  is  used  to  calculate  the  axial  position  of  the  sheath  during  the 
rundown  phase.  The  dynamic  inductance  of  the  sheath  is  dependent  on  the  positional 
tracking  of  the  arc. 

This  dynamic  inductance  expressed  as  Lpp  in  the  equivalent  circuit  is  dependent  on 
the  axial  position  z  as  well  as  the  distance  between  cathode  and  anode.  Utilizing  Faraday’s 
law  of  induction  for  a  single-turn  current  carrying  coil,  one  obtains  the  following 
expressions: 


d)  =  LI , 


where 


Be  z  dr , 


and 


B 


(31) 


(32) 

(33) 


Combining  the  previous  equations  will  give  a  general  solution  of  the  sheath  inductance  for 
the  device. 


L(t)  = 


ItpZ 
2  7t 


(34) 


This  general  inductance  is  adjusted  by  adding  an  exu-a  term  that  accounts  for  the  inductance 
due  to  the  radial  liftoff  of  the  plasma  sheath  from  the  insulator. 

The  time  dependent  z(t)  term  in  the  brackets  represents  the  inductance  over  the  anode,  while 
the  zio  term  is  the  contribution  of  the  sheath  inductance  over  the  insulator  region  during  the 
radial  liftoff  of  the  plasma  sheath.  The  term  z(t)  is  the  time  dependent  location  of  the 
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leading  edge  of  the  current  sheath  while  the  term  zio  is  the  length  of  the  insulator  and  thus  a 
fixed  quantity.  The  zio  contribution  to  the  inductance  describes  the  initial  inductance  before 
the  plasma  sheath  lifts  off  from  the  insulator  and  starts  to  propagate  down  the  device. 
From  this  total  expression  for  the  sheath  inductance,  the  initial  sheath  inductance  Lppft  =  0) 
can  be  found  as 


LpF(l  =  0)  =  ^ln(t)[A+z,o]. 

At  time  zero  the  time  dependent  expression  z(t)  can  be  replaced  by  the  sheath  thickness  A 
which  is  assumed  to  be  0.2  cm  in  the  transient  code. 

The  plasma  resistance,  Rpp.  is  dependent  on  the  temperature  as  well  as  the  density  of  the 
sheath.  The  resistance  is  derived  by  finding  the  electric  field  between  the  anode  and 
cathode  surfaces  using  Gauss's  law; 


I 


E  •dS  = 


£o 


with 


0^ 


Q 


27t  r  z  eo 


(37) 


(38) 


The  next  step  is  to  utilize  Ohm's  law  and  the  expression  for  the  current  through  the  plasma 
sheath 


R  =  ^, 
I 


where 


(39) 


(40) 


(41) 


Combining  these  three  equations  and  solving  for  the  resi.stance  give  the  final  expression  for 
resistance  as 
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(42) 


The  resistivity  t|  is  taken  to  be  the  Spitzer  resistivity  for  a  plasma  and  A  is  the  thickness  of 
the  plasma  sheath.  Plasma  is  considered  a  very  good  conductor  and  the  typical  values  of 
plasma  resistivity  are  very  small  in  magnitude  when  compared  to  external  circuit 
parameters. 

The  initial  plasma  sheath  resistance  was  calculated  using  the  initial  fill  gas  density 
and  an  assumed  initial  breakdown  temperature.  This  resistance  was  kept  constant  through 
the  simulation  and  has  negligible  effect  on  current  history. 

Before  the  plasma  sheath  begins  to  propagate  down  the  anode,  it  must  first  detach 
itself  from  the  surface  of  the  insulator.  Particle  pressures  resulting  from  the  arc  formation, 
anchor  the  sheath  to  the  insulator  until  the  applied  electromagnetic  force  overcomes  the 
ambient  pressure.  This  radial  lift-off  cannot  occur  unless  the  sheath  current  reaches  a 
particular  value.  One  can  describe  the  threshold  condition  as: 


2  P-o 


=  nkT. 


(43) 


this  is  the  balance  condition  between  the  magnetic  pressure  of  the  sheath  and  the  plasma 
pressure.  The  magnetic  induction  B  can  be  described  by  a  previous  expression  derived 
from  Ampere’s  law  for  a  straight  current  carrying  conductor  as  in  Equation  (33).  Equation 
(33)  is  substituted  into  the  pressure  balance  (Equation  43),  and  solved  for  the  current  I. 
The  result  is  the  expression  for  calculating  liftoff  current. 

I  =  ^/IZSISkT 

V  Po  A  ’  (44) 

Where  kT  is  the  plasma  sheadi  temperature  in  joules,  Na  is  Avogadro's  number,  p  is  the 
initial  fill  gas  density,  A  is  the  atomic  mass  of  the  fill  gas  and  rg  is  the  anode  radius. 

In  order  to  calculate  a  lift-off  current,  an  initial  sheath  temperature  must  be  assumed 
for  determining  the  particle  pressure.  It  is  assumed  that  the  sheath  pressure  remains 
constant  during  the  rundown  phase.  Significant  heating  of  the  sheath  is  assumed  to  occur 
during  the  radial  collapse  of  the  sheath  due  to  compressional  heating,  and  not  in  the 
rundown  phase.  This  lift-off  current  is  denoted  as  Ilo  ^4  is  responsible  for  a  delay  in  the 
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initiation  of  plasma  rundown.  In  order  to  account  for  this  lift-off  delay,  the  lift-off  current 
was  added  to  the  Snowplow  model  previously  derived.  The  resultant  expression  gives 


z 


z 


(45) 


for  the  modified  snowplow  velocity. 

The  previously  derived  plasma  relations  are  directly  coupled  to  the  circuit 
parameters  in  the  load  branch  representing  the  plasma  sheath.  These  sets  of  coupled 
equations  realistically  model  the  current  history  and  rundown  velocity  of  the  dense  plasma 
focus  device. 


Simulation  of  the  Livermore-I  Experiment 

Once  the  equations  describing  the  plasma  sheath  properties  were  developed  and 
implemented  into  the  transient  solver,  the  code  was  used  to  simulate  the  Livermore-I 
plasma  focus  experiment.  The  input  parameters  and  electrode  geometry  of  the  Livermore-I 
device  were  inserted  into  the  input  deck  of  the  code.  Table  2  shows  the  input  parameters 
for  this  validation  test.  The  calculated  results  from  the  code  were  plotted  and  when 
possible,  compared  against  the  experimental  values.  The  Livermore-I  experiment  measured 
the  total  external  current  of  the  equivalent  circuit  and  an  inferred  leakage  current.  This 
leakage  current  has  not  been  directly  measured,  but  it  is  indicative  of  a  current  loss 
mechanism. 

Utilizing  the  anode  profile  of  the  Livermore-I  experiment  (Figur''  15),  Figure  16 
shows  the  results  of  the  calculated  external  and  leakage  currents  and  contrasts  them  with  the 
experimental  results  from  Livermore-I  experiment.  There  is  good  agreement  between  the 
calculated  and  experimental  curves,  although  the  code  underpredicts  external  current  at  the 
start  of  the  transient  and  overpredicts  current  at  the  end  of  the  run.  This  is  thought  to  be 
due  to  the  distortion  of  the  sheath  at  lift-off  and  at  collapse.  The  profile  of  the  plasma 
sheath  would  determine  the  flux  area  between  conducting  surfaces  and  thus  the  sheath 
inductance.  Variations  in  the  inductance  would  affect  the  current  history  of  the  sheath,  and 
it  would  appear  that  assuming  a  sheath  profile  perpendicular  to  the  anode  is  good  for  times 
between  the  lift-off  phase  and  when  the  sheath  collapses  to  the  axis. 
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Table  2.  Input  Parameters  for  the  Livermore-I  Plasma  Focus  Experiment 


Input  Parameters 

Livermore-I  Experiment 

Voltage  (volts) 

2.7  X  104 

Initial  Inductance  (henries) 

2.5  X  10-8 

External  Resistance  (ohms) 

5.0  X  10-3 

Leak  Resistance  (ohms) 

1.2  X  10-1 

External  Capacitance  (farads) 

3.55  X  10-4 

Anode  Radius  (cm) 

5.08 

Cathode  Radius  (cm) 

8.0 

Anode  Tip  Radius  (cm) 

1.27 

Length  of  Insulator  zio  (cm) 

14.0 

Initial  Point  of  Anode  Curvature  (cm) 

31.4 

■Tip  of  Anode  z^p  (cm) 

38.2 

Fill  Gas  Density  (g/cm^) 

2.2  X  10-7 

Atomic  Weight  of  Fill  Gas  (amu) 

1 

Assumed  Lift-off  Temperature  (eV) 

5 

Time  Step  (sec) 

5.0  X  10-8 
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Anode  Radius  (m)  {  x  10 ' 


Figure  15 

The  Anode  Profile  of  the  Livermore-I  Device 
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Figure  16 

Experimental  and  Calculated  Current 
Histories  for  the  Livermore-I  Experiment 


Figures  17  through  19  show  the  calculated  values  of  plasma  sheath  inductance,  axial 
rundown  velocity,  and  the  node  voltage  across  the  sheath  and  leakage  components. 
Experimental  results  were  not  available  for  these  plots  due  to  the  difficulty  involved  in 
measuring  these  quantities[21].  However,  the  calculated  results  represented  the  proper 
trend  in  behavior  for  the  rundown  phase  of  operation.  The  axial  velocity  in  Figure  18 
remains  at  zero  until  the  sheath  current  reaches  the  lift-off  value  and  begins  to  propagate 
down  the  anode.  Axial  velocity  also  displays  a  slowing  down  trend  as  the  sheath  hits  the 
curved  portion  of  the  anode(see  Figure  15),  then  it  experiences  an  acceleration  as  the  anode 
straightens  out.  This  matches  correctly  with  the  inductance  behavior  at  lift-off  and  at  the 
end  of  rundown.  The  inductance  remains  constant  before  the  sheath  lifts  off  of  the 
insulator  and  increases  as  the  cross  sectional  flux  area  between  the  conductors  increases. 
The  voltage  profile  in  Figure  19  also  shows  the  voltage  behavior  expected  during  the 
transient.  In  the  early  stages  of  the  transient,  the  gap  voltage  ramps  up  due  to  the  effect  of 
the  fast  rising  current  on  the  constant  plasma  sheath  inductance.  The  voltage  plateau  in  the 
middle  of  the  transient  is  due  to  the  rising  inductance  that  occurs  after  the  plasma  sheath  has 
started  to  propagate  axially  down  the  anode.  At  the  end  of  the  transient,  the  inductance 
increase  because  of  the  increase  in  cross  sectional  flux  area  that  is  due  to  the  tapering  anode 


Figure  17 

Calculated  Plasma  Sheath  Inductance  for  the 
Livermore-I  Experiment 
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Time  (sec)  (  x  10^  ) 
Figure  18 

Calculated  Axial  Rundown  Velocity  for 
the  Livennore-I  Experiment 


tip.  This  last  increase  in  inductance  accounts  for  the  voltage  increase  that  occurs  at  the  end 
of  the  transient. 

The  calculated  results  from  the  transient  code  were  in  good  agreement  with  the 
experimental  results  and  the  general  behavior  that  was  expected.  The  differences  in  current 
history  at  the  beginning  and  ending  of  the  transient  were  due  to  the  radial  behavior  of  the 
sheath  at  liftoff  and  at  the  onset  of  sheath  collapse.  These  are  two-dimensional  effects  that 
cannot  be  easily  accounted  for  using  the  one-dimensional  snowplow  model.  The 
inductance,  node  voltage,  and  axial  velocity  behaved  as  expected  and  it  should  be  noted 
that  the  snowplow  model  predicted  a  rundown  velocity  that  did  not  exceed  the  implosion 
velocity  limit  of  3.5  x  10^  m/s.  This  limit  is  based  on  the  implosion  velocity  of  an  inertial 
confinement  fusion  target  and  is  a  good  physical  limit  for  the  maximum  velocity  attainable 
for  the  plasma  sheath. 
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Figure  19 

Calculat'd  Gap  Voltage  Across  the  Plasma  Sheath 
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PARAMETRIC  STUDY  OF  TOE  DPF  ELECTRODES 


Introduction 

Parametric  tcjjting  of  the  effects  of  electrode  configuration  is  addressed  in  this 
section.  Initial  studies  on  the  application  of  the  dense  plasma  focus  as  a  .space  propulsion 
concept  were  done  on  the  assumption  that  a  device  similar  to  the  Livcrmore-I  experiment 
could  be  made  to  produce  the  necessary  current  for  fusion  ignition.  The  purpose  of  this 
parametric  study  is  to  make  a  realistic  as.sessment  on  this  assumption  that  was  made 
previously.  In  order  to  accomplish  this  objective,  a  variation  of  the  radial  dimensions  of 
the  anode  will  be  implemented.  All  of  the  initial  input  parameters  other  than  the  radial 
dimensions  will  be  kept  constant  and  the  effects  of  these  electrode  changes  will  be 
documented.  If  the  radial  variation  does  not  produce  the  required  ignition  current, 
additional  radial  and  axial  vanations  will  be  implemented.  Once  the  required  ignition 
current  is  reached,  the  results  from  the  end  of  rundown  will  be  input  into  a  code  which 
calculates  the  performance  of  a  DPF  propulsion  system.  The  system  performance  will  be 
contrasted  to  previous  system  calculations  that  were  carried  out  with  assumed  rundown 
values.  Again,  the  goal  of  this  work  is  to  provide  a  realistic  as.sessment  of  the  parameter 
requirements  that  are  needed  to  make  the  DPF  device  feasible  as  a  space  propulsion  system. 

Radial  Variation  of  Anode 

Plasma  focus  sheath  current  for  the  Livermore-I  type  device  is  sensitive  to  changes 
in  the  annular  gap  distance  between  cathode  and  anode.  Reduction  of  this  gap  will  lead  to  a 
reduced  inductance  which  increases  the  magnitude  of  the  plasma  sheath  current.  The  gap 
reduction  will  also  reduce  the  volume  of  fill  gas  contained  in  the  annulus.  This  means  that 
less  gas  will  be  entrained  in  the  sheath  as  it  propagates  down  the  device,  which  will  limit  the 
density  in  the  pinch  region.  Ideally,  the  current  and  density  at  the  end  of  the  rundown  phase 
should  be  optimized  to  present  the  best  po.ssible  conditions  prior  to  the  pinch  phase. 

The  radial  parameter  tests  were  conducted  by  increasing  the  anode  radius  in  order  to 
get  a  larger  sheath  current.  The  gap  distance  between  the  cathode  and  anode  were  reduced 
by  half  for  each  test  case  and  the  resulting  currents  were  documented  in  Table  3. 
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Four  different  radial  variations  of  the  anode  were  used  in  these  tests.  The  gap 
distances  were  reduced  from  the  initial  Li’.'ermore-I  value  of  0.0292  meter  down  to 
0.00365  meter.  Figure  20  shows  the  different  profiles  of  the  anode  for  these  variations. 
Anode  length  and  other  plasma  focus  parameters  were  kept  identical  to  the  initial 
Livermore-I  configuration.  Figures  21  through  24  show  the  sheath  currents,  inductances, 
electrode  gap  voltages  and  rundown  velocities  of  each  test  set.  The  maximum  sheath 
current  is  reached  at  time  tmax  •  It  can  be  seen  from  the  curves  in  Figure  2 1 ,  that  the 
maximum  sheath  current  corresponding  to  the  smallest  gap  length  is  only  1.527  MA.  This 
is  far  below  the  required  current  needed  to  produce  an  adequate  pinch  temperature.  Further 
reduction  of  the  gap  distance  would  result  in  an  unacceptably  low  numbc  j.  nsity  in  the 
pinch. 

Despite  changes  in  the  anode  geometry,  it  is  apparent  that  the  DPF  in  the 
Livermore-I  configuration  cannot  supply  the  magnitude  of  current  needed  for  fusion 
ignition  due  to  the  large  electrode  gap.  Reduction  of  the  electrode  gap  distance  will  increase 
sheath  current  but  also  decrease  total  annular  volume.  The  reduction  in  annular  volume  will 
mean  that  fewer  fuel  gas  particles  v'ill  be  entrained  into  the  pinch  region  after  rundown. 
We  can  conclude  that  a  Livermore-I  type  device  is  unsuited  to  the  task  of  producing 
currents  of  10  or  20  MA.  This  conclusion  contrasts  sharply  to  the  previous  sumption 
that  a  Livermore-I  type  device  could  produce  this  large  magnitude  of  current 

Variations  of  Elecuode  Length  and  Charging  Voltage 

In  order  to  facilitate  the  attainineat  of  a  10-20  MA  sheath  current,  device  parameters 
other  than  the  anode  radius  must  be  adjusted.  The  charging  voltage  of  the  capacitor  bank 
will  be  increased  in  concert  with  the  differing  anode  radii.  The  increases  in  charging 
voltage  are  conducted  with  the  assumption  that  the  initial  external  inductance  and  resistance 
remain  at  the  same  values.  Lack  of  experimental  information[22]  on  such  high  current 
circuits  necessitates  the  need  for  this  assumption.  The  initial  charging  voltage  is  increased 
for  each  test  in  increments  of  27  kV,  each  new  voltage  is  tested  over  three  differing  radial 
lengthsfO. 0508m,  0.0654m,  and  0.0727m).  Boosting  the  voltage  increases  the  current  in 
the  sheath  and  thus  the  rundown  velocity.  This  increase  in  axial  velocity  reduces  the 
rundown  time  of  the  sheath.  As  a  result,  the  sheath  runs  out  of  the  device  before  peak 
current  can  be  reached.  In  order  to  counteract  this,  the  electrode  length  is  accordingly 
increased  so  that  the  sheath  runs  out  of  the  device  when  the  peak  current  in  the  plasma  arc 
is  reached. 
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Sizing  of  the  electrode  length  maximizes  the  geometry  of  the  device  for  each  new  voltage, 
radius  combination. 

The  voltage,  and  electrode  variational  tests  were  conducted  until  cune.-iis  of  10, 1 
and  20  MA  were  obtained  for  each  separate  test  combination.  Table  3  shows  the  optimized 
voltage  and  geometry  needed  to  obtain  these  necesssary  currents. 


Table  3.  Optimized  Parameters  for  High  Current  Delivery 


Electrode  Gap 

Anode 

Cathode 

Electrode 

Maximum 

Distance  (m) 

Radius  (m) 

Length  (m) 

Current  (MA) 

mmKsmm 

0.08 

■■RHiaHl 

0.08 

243 

bbi^h 

0.08 

243 

0.0727 

0.08 

0.08 

378 

0.08 

1.952 

15.04 

*324 

0.0727 

1.232 

BB&3B 

378 

0.0654 

1.482 

486 

0.0508 

0.08 

3.102 

BBESB 

*  Best  Case 

It  is  desirable  to  minimize  the  voltage  requirements  as  well  as  the  geometrical 
dimensions  of  this  device.  For  this  reason,.the  minimum  voltage  cases  for  the  10,15,  and 
20  MA  range  shall  be  used  as  the  final  voltage,  geometry  combinations.  An  assumption  is 
made  that  the  arc  current  will  not  saturate  in  the  inter-electrode  gap  during  operation. 
Another  assumption  made  is  that  the  electrode  material  can  withstand  the  high  temperature 
created  by  joule  heating  and  the  plasma  arc.  Each  of  the  final  chosen  geometries  have  gap 
widths  of  0.0073  meter,  this  allows  the  anode  radius  to  be  increased  from  the  previous 
Livermore-I  geometry.  The  new  anode  radius  gives  an  increased  electrode-arc  interface 
area  which  results  in  a  decreased  current  density.  The  decrease  in  effective  current  density 
is  an  added  advantage  for  the  chosen  final  geometries.  The  best  case  voltage-electrode 
combination  is  marked  by  an  asterik  in  Table  3.  This  combination  was  chosen  because  it 
utilized  the  lowest  charging  voltage  for  the  20  MA  cases.  However,  even  the  best  case 
scenario  utilizes  an  extremely  high  magnitude  of  input  voltage  which  would  require  a  Marx 
generator  configuration  for  the  input  circuit.  The  20  MA  case  is  considered  to  be  optimal 
for  the  DPF  in  the  space  propulsion  application  since  the  high  current  will  provide  higher 
plasma  temperatures  and  enhanced  reaction  rates. 
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DENSE  PLASMA  FOCUS  PROPULSION  SYSTEM 


Introduction 

The  use  of  the  dense  plasma  focus  as  a  viable  propulsion  concept  requires  that  the 
system  design  not  exceed  practical  standards  for  operation.  The  parametric  study 
conducted  in  this  work  is  a  facet  of  DPF  system  design  that  was  not  rigorously  explored 
previously  by  Choi  and  Leakeas[  1].  The  transient  code  calculations  obtained  in  this  work 
will  be  fed  into  the  DPF  system  code  [1].  The  basic  system  design  and  requirements  will 
be  kept  identical  to  the  previously  tested  model,  the  only  exception  being  the  elimination  of 
the  assumption  that  a  Livermore  I  type  device  could  be  used  to  generate  the  necessary 
current.  The  system  code[l]  will  be  rerun  for  the  modified  electrode  geometries  and 
charging  voltages  that  were  obtained  from  the  parametric  study  in  the  previous  chapter. 
The  following  subsections  contain  basic  descriptions  of  the  guiding  principles  used  in  the 
design  of  the  DPF  system  code. 

Rocket  Dynamics 

The  performance  characteristics  of  a  rocket  propulsion  system  can  be  judged  using 
several  important  parameters  which  describe  the  power  and  efficiency  of  a  system.  One  of 
the  most  vital  parameters  that  describes  rocket  performance  is  specific  impulse  or  Isp. 
Specific  impulse  is  defined  as  the  amount  of  momentum  gained  per  weight  of  fuel  burned, 
and  is  given  by  the  expression 


I 


(46) 


where  Vex  is  the  exhaust  velocity,  and  g  is  the  gravitational  constant  of  earth.  The  specific 
impulse  uses  the  gravitational  constant  of  earth  instead  of  the  local  gravitation  because  it  is  a 
unif  d  reference  value  for  all  types  of  propulsion  devices.  If  one  were  to  use  the  local 
gravitational  constant  in  deep  space,  the  gravity  would  be  very  minute  and  result  in 
unreasonably  high  values  of  specific  impulse.  Therefore,  it  is  expedient  to  normalize  the 
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weight  of  fuel  to  the  gravitational  constant  of  the  point  of  origin.  Specific  impulse  is 
measured  in  units  of  seconds  and  is  a  good  measure  of  the  efficiency  of  a  rocket.  Thrust  is 
another  important  performance  parameter  of  a  rocket  system.  The  expression  for  thrust  is 
given  by  the  following  expression 


^thrust  ~  ^propellant  '^ex  ^Pex  "  Pamb^ 


(47) 


where  F  is  expressed  in  units  of  Newtons.  The  second  term  on  the  right  of  the  equation 
represents  the  pressure  force  due  to  the  differential  pressure  between  the  exhaust  stream 
and  the  ambient  pressure.  This  term  is  taken  to  be  negligible  in  comparison  to  the  primary 
force  contributed  by  the  rocket  exhaust  flow.  The  exhaust  power  of  a  rocket  can  be 
expressed  in  terms  of  thrust  and  specific  impulse  and  is  given  by 


P  =  ^  F  V, 
2 


ex  ~  2  S  P  ^sp 


(48) 


This  relation  implies  that  for  a  fixed  power,  any  increase  in  specific  impulse  will  deman<^  a 
similar  decrease  in  thrust.  In  order  to  optimize  a  propulsion  system,  both  of  these 
parameters  must  be  maximized.  The  last  parameter  necessary  for  defining  performance  is 
the  burnout  velocity  of  the  vehicle  or  Av.  The  Av  is  derived  from  the  equation  of  motion  in 
free  space  for  a  rocket  and  is  given  by 


Av  =  v„  In^, 


m 


(49) 


where  Av  is  the  velocity  increment,  mo  is  the  initial  mass  of  the  entire  vehicle  including  fuel 
and  payload,  and  m  is  the  final  mass  of  the  vehicle  without  fuel.  These  parameters  will 
play  a  key  role  in  the  understanding  of  the  effectiveness  of  the  DPF  propulsion  system. 


Fusion  Principles 

The  Dense  Plasma  Focus  thruster  relies  on  the  generation  of  fusion  power  to 
provide  vehicle  thrust.  Fusion  is  a  thermonuclear  reaction  involving  the  fusing  of  ions,  and 
is  capable  of  generating  large  amounts  of  energy  per  unit  volume.  In  order  for  fusion  to 
occur,  the  ions  must  be  forced  together  in  order  to  overcome  the  Coulomb  repulsive  force 
that  naturally  repels  particles  with  like  charges.  To  accomplish  this,  a  collection  of  high 
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temperature  charged  particles,  also  known  as  a  plasma,  must  be  confined,  compressed  and 
heated. 

Fusion  reactions  depend  largely  on  the  plasma  temperature,  density,  and 
confinement  time.  An  expression  for  fusion  power  density  is 


pF  =  n,  n2(cr  v)Wo  ^  (50) 

where  ni  and  n2  are  the  number  densities  of  the  non-equal  reacting  species,  is  the 
reaction  rate  of  the  plasma,  and  Wq  is  the  energy  liberated  per  reaction.  The  reaction  rate 
parameter  is  a  temperature  dependent  quantity  which  measures  how  quickly  a  reaction  takes 
place  at  a  given  temperature.  In  order  to  calculate  the  total  energy  yield  from  a  constant 
volume,  simply  multiply  the  power  density  by  the  volume  of  plasma  and  by  the  time  over 
which  the  reaction  extends. 

Since  these  parameters  are  dependent  on  the  choice  of  fusion  fuel  used,  it  is 
advantageous  to  select  an  optimal  fuel  for  a  specific  application.  For  the  case  of  the  DPF 
thruster,  we  would  like  a  fuel  that  fulfills  certain  criteria; 

1 .  Provides  a  high  energy  output  per  reaction. 

2.  Maximizes  the  reaction  rate  parameter  at  "low"  operating  temperatures  (keV  range). 

3.  Suppresses  neutron  production  since  neutrons  cannot  be  directed  ^ith  a  magnetic 
field. 

Several  fuel  choices  were  considered  as  possible  candidates  for  DPF  fuel,  these  fuels  are 
listed  below. 


D*  +  V-^  4He^+  (3.5MeV)  +  n  (14.1MeV)  ^ 

(l.OlMeV)  +  p  (3.02MeV)  (50%) , 

(.82MeV)  +  n  (2.45MeV).  (50%) , 

D+  +  He^-^  ^He-^  (3.6MeV)  +  p  (14.7MeV)  . 


(51) 

(52) 

(53) 


The  first  reaction  listed  is  a  deuterium-tritium  reaction,  second  is  the  DDn  and  DDp 
reaction,  and  lastly  the  deuterium-helium-3  reaction.  The  D-T  reaction  has  the  highest 
reaction  rate  at  low  temperature  among  three  reactions,  but  80%  of  the  reaction  energy  is 
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carried  away  by  the  neutron.  This  is  considered  detrimental  to  the  production  of  thrust, 
since  the  neutrons  will  fly  in  any  direction  and  cannot  be  channeled  by  a  magnetic  field. 
The  D-D  reaction  is  split  into  two  sub-reactions(DDn,DDp),  each  reaction  having  an  equal 
probability  of  occurring.  This  reaction  is  also  not  desirable  since  there  is  a  .50%  probability 
of  producing  a  neutron.  This  leaves  D-^He  as  the  remaining  candidate  for  the  "low ' 
temperature  fuels.  The  D-^He  reaction  produces  the  highest  amount  of  energy  per  reaction 
(18.3MeV)  and  does  not  produce  neutrons  in  its  primary  reaction.  The  reaction  rate  for  a 
D-^He  reaction  is  also  comparable  to  the  D-T  reaction  rate  at  similar  temperatures  (keV 
range).  Although  D-^He  does  not  produce  any  neutrons,  secondary  neutron  production  in 
a  D-^He  fuel  is  possible  from  background  D-D  reactions.  The  initial  DPF  system  study  by 
Choi  and  Leakeas  [1],  considered  other  possible  advanced  fuels  such  as  proton-Lithium-6 
and  proton-Boron-1 1,  but  these  fuels  require  ignition  temperatures  beyond  reasonable 
limits  (past  100  keV)  and  were  discarded  as  possible  choices. 


Brief  Description  of  DPF  Propulsion  System  Code 

The  DPF  system  code  (not  to  be  confused  with  the  transient  code  developed  in  this 
work)  calculates  the  output  power  of  the  pinch  and  the  resulting  system  performance  of  a 
DPF  propulsion  system.  Many  input  parameters  were  assumed  during  the  initial  DPF 
study  conducted  previously[l].  The  rundown  calculations  in  the  system  code  were 
discarded  and  replaced  by  the  circuit  transient  code.  This  was  done  in  order  to  provide  a 
more  realistic  assessment  of  the  design  requirements  of  the  DPF.  However,  the  pinch 
phase  calculations  and  system  calculations  are  unchanged  from  the  previous  system  code. 

The  following  model  provided  a  basis  for  the  parameter  calculations  in  the  pinch 
region.  The  pinch  was  modeled  as  a  cylindrical  region  of  assumed  radius  and  length  that 
contained  a  certain  fraction  of  sheath  plasma.  In  order  to  determine  the  temperature  inside 
the  pinch  region,  a  balance  between  magnetic  pressure  and  plasma  pressure  was  assumed. 


npkT  =  :^ 
2  Po 


(54) 


The  azimuthal  magnetic  induction  at  the  surface  of  the  pinch  volume  is  Bq,  Po  is  the 
permeability  of  free  space,  np  is  the  pinch  number  density  and  kT  is  the  product  of 
Boltzmann's  constant  and  the  plasma  temperature  in  degrees.  The  pinch  number  density  np 
was  assumed  to  be  a  fraction(f)  of  the  initial  gas  density  present  in  the  annular  region.  This 
trapping  fraction  was  chosen  in  order  to  provide  a  match  with  experimental  values  from 
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Livermorc-I  data  (np  -  Expressing  np  in  terms  of  the  electrode  dimensions  and 

initial  fill  gas  density  gives 


_  f  Pi  la(rc  -  rj) 

"P  ~  I  7 

*P  ‘'p  "^P  (55) 

where  la  and  Ip  are  the  anode  and  pinch  lengths,  respectively,  rq,  ra  and  rp  are  the  radius  of 
the  cathode,  anode  and  pinch,  and  mp  is  the  average  mass  of  particles  in  the  pinch. 
Utilizing  Ampere's  Law  and  integrating  around  the  cylindrical  surface  of  the  pinch  region 
gives  an  expression  for  the  azimuthal  magnetic  induction  Be  (as  equation  33) 

2  "  '■p  (56) 

By  combining  Equations  (  54)  through  (  56 ),  a  final  equation  relating  plasma  temperature 
to  current  is  obtained. 


IcT  —  ^P  ^P 

8  f  pi  U  (ri  -  ri)  (57) 

The  pinch  analysis  model  assumes  that  this  scaling  holds  regardless  of  input  current 
magnitude. 

The  system  code  relied  on  the  pinch  parameters  as  an  input  for  calculating  the 
resultant  fusion  power.  All  of  the  performance  parameters  of  the  system  are  calculated 
around  this  energy  yield.  The  DPF  system  itself  consists  of  the  dense  plasma  focus  device, 
storage  tanks  for  fuel  and  coolant,  capacitor  and  charging  circuit  mechanisms  and  a 
magnetic  nozzle.  This  standard  configuration  of  the  DPF  system  is  displayed  in  Figure  25. 
The  hydrogen  tank  supplies  the  necessary  cryogenic  needed  to  cool  the  walls  of  the  cathode 
and  the  combustion  chamber.  This  hydrogen  can  also  be  injected  into  the  exhaust  flow  of 
the  rocket  to  increase  the  resultant  thrust  The  capacitor  banks  are  used  to  provide  the 
discharge  current  necessary  for  operation  and  the  turbine-generator  is  used  to  re-energize 
the  banks  after  each  discharge  pulse.  The  turbine  is  driven  by  the  cryogenic  coolant  after  it 
has  cooled  the  electrode  and  combustion  chamber  walls.  The  magnetic  nozzle  channels  the 
exhaust  flow  of  particles  out  of  the  combustion  chamber[21].  The  primary  purpose  of  the 
magnetic  nozzle  is  to  prevent  the  fusion  products  from  impacting 
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Figure  25 

System  Diagram  for  Dense  Plasma  Focus  Thruster 


with  the  wall  material  as  well  as  generating  and  accelerated  exhaust  tlow.  For  the 
continuous  and  impulsive  modes  of  firing,  the  fusion  products  are  diluted  with  the  added 
hydrogen  intlow  which  attenuates  the  charged  panicle  temperature. 

However,  it  is  still  desurable  to  keep  the  high  energy  flow  away  from  the  wails  which  will 
degrade  the  structural  integrity  of  the  combustion  chamber. 

An  important  factor  forjudging  the  performance  of  the  DPF  propulsion  system  is 
the  total  system  mass.  If  this  mass  is  too  high,  thrust-weight  ratios  are  reduced,  thereby 
cutting  performance.  Capacitor  mass  is  a  dominant  component  of  the  overall  system  ma.s.s 
and  any  reduction  of  this  is  highly  desirable.  Current  technology  allows  a  specific  energy 
of  about  0.2  kJ/kg  for  modem  capacitors.  In  order  to  fulfil!  the  requirements  for  the 
previous  system  model,  capacitor  masses  on  the  order  of  40,{KX)  kg  would  be  needed.  To 
offset  this  problem,  a  further  assumption  was  made  that  capacitors  with  specific  energic.s  of 
2  kJ/kg  could  be  obtained!  1].  This  would  reduce  the  capacitor  ma.ss  considerably  and 
allow  greater  thrust-to-weight  ratios.  The  system  capacitor  mas.se.s  were  previously 
calculated  using  an  l2  scaling  law  which  relates  the  ratio  of  capacitor  mas.ses  to  a  ratio  of 
squared  currents.  The  scaling  law  is  derived  by  equating  the  ratio  of  capacitor  masses  to 
the  ratio  of  energy  expended  in  the  operation  of  the  device 

M  -  W 

M„  "  Wo  ■  (58) 


The  terms  Mq  and  Wq  are  the  capacitor  mass  and  expended  energy  of  a  base  case 
experiment  while  M  and  W  are  the  new  device  values.  The  expended  energy  in  the  den.se 
plasma  focus  can  be  expressed  as  the  product  of  the  magnetic  force  dnving  the  sheath  and 
the  sheath  displacement 


2  iio  '  ’  (59) 

p  —  ^ 

“e  -  y  'z:_  ’ 
where  ^  ^ 

and  Vp  is  the  radius  of  the  pinch  region  Substituting  Be  and  Equation  (59)  into  (58)  yields 
the  scaling  expression  for  capacitor  mass. 
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lo  (r?  -  tIo)  Ay 


(60) 


The  previous  DPF  system  study  [1]  assumed  that  the  Livermore-I  base  case  configuration 
would  not  have  to  be  changed  in  order  to  attain  the  current  necessary  for  operation.  This 
assumption  eliminated  the  geometrical  dependence  in  the  scaling  law  since  both  the  base 
case  and  new  configuration  had  identical  dimensions.  However,  this  study  will  utilize  the 
geometrical  dependence  in  the  scaling  law  since  the  test  case  geometries  differ  from  the 
Livermore-I  base  case  geometry. 

Another  important  factor  that  influences  performance  is  the  mode  of  operation  of 
the  system.  Choi  and  Leakeas  [1]  considered  three  different  modes  of  operation:  I  ) 
operation  of  the  DPF  as  a  closed  system  with  no  addition  of  hydrogen  into  the  exhaust 
stream,  2)  continuous  firing  with  addition  of  hydrogen,  3)  firing  for  short  periods  of  time 
with  large  exhaust  of  hydrogen.  The  third  mode  was  found  to  be  the  most  advantageous 
since  the  large  exhausts  of  hydrogen  increased  the  exhaust  mass  flow  rate  and  maximized 
thrust-to- weight  ratio. 

Performance  Results  with  Modified  Electrode  Configurations 

The  new  electrode  configurations  developed  in  Chapter  3  were  tested  on  the  DPF 
system  code  and  contrasted  with  the  results  from  the  previous  system  study  [1],  The  tests 
were  run  for  differing  values  of  Av  using  the  enhanced  electrode  configurations.  Plots 
were  obtained  for  the  specific  impulse  (Isp)  and  thrust-to-weight  ratio  for  each  Av 
requirement.  High  Av  requirements  are  necesssary  in  order  to  shorten  trip  time  for  longer 
range  missions.  The  reduced  trip  time  will  minimize  the  vehicle  occupants  exposure  to  zero 
gravity  and  cosmic  radiation.  The  mode  of  operation  used  in  these  test  scenarios  is  the 
impulsive  firing  mode.  It  is  one  of  the  three  operational  modes  previously  mentioned,  and 
it  involves  pulsing  the  thruster  while  exhausting  large  amounts  of  hydrogen  into  the 
exhaust  flow.  The  advantages  of  this  mode  of  operation  is  that  it  greatly  reduces  system 
mass  due  to  the  exhausting  of  massive  amounts  of  hydrogen,  this  in  turn  increases  the 
thrust-to-weight  ratio  significantly.  The  increased  thrust-to-weight  ratio  allows  the  vehicle 
to  accelerate  rapidly  until  the  fuel  is  exhausted. 

The  enhanced  electrode  configurations  take  advantage  of  the  capacitor  mass  scaling 
law  shown  previously.  The  I^  scaling  with  geometrical  dependence  predicts  smaller 
capacitor  masses  than  the  previous  P  scaling  without  geometrical  dependence.  The  result 
is  that  the  capacitor  mass  needed  for  the  enhanced  configurations  is  smaller  than  that  for  a 
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Livermore-I  type  geometry.  A  nominal  base  case  was  chosen  in  the  previous  DPF  system 
study.  This  base  case  was  for  a  plasma  sheath  operating  current  of  20  MA  and  a  Av 
requirement  of  10  km/s.  Table  4  shows  the  comparison  between  the  base  case  outputs  for 
the  assumed  Livermore-I  coni'iguration  and  the  enhanced  20MA  elecu'ode  configuration  of 
this  study. 

Figures  26  through  3 1  show  the  plots  for  specific  impulse  and  thrust-to-weight 
ratios  for  different  Av  requirements.  For  Av=5km/s  (Figure  26),  the  specific  impulse  of  the 
Livermore-I  geometry  electrodes  and  the  enhanced  electrodes  showed  little  difference.  The 
best  case  appears  to  be  for  the  20MA  enhanced  electrode  configuration.  The  20MA  current 
provided  a  better  fusion  bum  in  the  pinch  which  resulted  in  a  higher  exhaust  velocity  over 
the  varying  range  of  propellant  mass  flow  rates.  The  propellant  mass  flow  rate  is  the  mass 
flow  rate  of  hydrogen  that  is  injected  into  the  combustion  chamber.  This  injection  of 
hydrogen  will  increase  the  thrust  of  the  device  by  increasing  the  mass  expelled  from  the 
exhaust  of  the  vehicle.  However,  if  the  mass  flow  rate  increases  for  a  fixed  current,  the 
exhaust  velocity  decreases  due  to  the  collisional  transfer  of  energy  between  the  the  charged 
particle  fusion  products  and  the  injected  hydrogen.  The  resultant  uend  is  the  degradation  of 
specific  impulse  for  increasing  hydrogen  mass  flow  rates.  Figure  27  shows  the  thrust-to- 
weight  ratio  for  the  Av=5km/s  case.  The  thrust-to-weight  ratio  of  the  enhanced 
configuration  shows  a  marked  improvement  over  the  previously  assumed  Livermore  I 
geometry.  This  plot  shows  the  trend  of  thrust-to-weight  ratio  increasing  as  the  propellant 
mass  flow  increases.  Again,  the  best  case  appears  to  occur  for  the  20MA  enhanced 
electrode  configuration. 

The  other  test  cases  for  Av=20  and  40km/s  show  basically  the  same  trend  as  the 
Av=5km/s  case.  The  enhanced  electrode  configuration  operating  at  20MA  seems  lo  provide 
the  optimal  case  for  each  Av  requirement.  For  increasing  Av  requirements  the  thrust-to- 
weight  ratio  curve  tends  to  flatten  out  or  decrease  for  increasing  propellant  mass  flow  rates. 
This  is  due  to  the  extra  propellant  mass  that  must  be  carried  in  order  to  reach  the  required 
Av.  The  increase  in  propellant  mass  will  decrease  the  thrust-to-weight  ratio  attainable  by 
the  system.  Inevitably,  we  come  to  the  same  conclusion  as  the  previous  study,  that  the 
DPF  system  thruster  is  most  efficient  for  low  Av  requirement  applications. 


50 


Table  4.  Propulsion  Parameters  for  Base  Case 
(Av=lOkm/s,  1=20MA) 


Livermore  I 

Enhanced 

Electrodes 

Electrodes 

Rundown  Velocity 

Vrun 

(mys) 

3.35x106 

3.5x105 

Deuterium  Bumup  Fraction 

fo 

0.699 

0.694 

Helium-3  Bumup  Fraction 

I'He 

0.442 

0.433 

D-^He  Fusion  Power 

PoHe 

(MW) 

2743.854 

2910.19 

DDn  Fusion  Power 

PDDn 

(MW) 

18.65955 

20.41 

DDp  Fusion  Power 

PODp 

(MW) 

83.55009 

91.79 

Totd  Fusion  Power 

Pf 

(MW) 

2846.066 

3022.40 

Power  to  Focus 

Pin 

(MW) 

5.4 

64.80 

Bremsstrahlung  Loss 

Pb 

(MW) 

15.79 

26.75 

Cyclotron  Loss 

Pc 

(MW) 

81.49 

118.6098 

Total  Power  Loss 

Pl 

(MW) 

97.2943 

145.3646 

Power  Increase 

AP 

(MW) 

2743.363 

2812.235 

Total  Mass  Flow 

Mt 

(kg/s) 

31.00 

31.233 

Propellant  Thrust 

Fp 

(N) 

4.73x105 

4.84x105 

Totk  Bum  Time 

tb 

(s) 

3479.949 

3328.988 

Payload  Mass 

Ml 

(kg) 

1x105 

lxl05 

Propellant  Mass 

Mp 

(kg) 

2.12x105 

2.046x105 

Propulsion  System  Mass 

Msys 

(kg) 

3.18x104 

3.06x104 

Fuel  Mass 

Mp 

(kg) 

0.7018 

0.59705 

Fuel  System  Mass 

Mpsys 

(kg) 

0.0702 

0.0597 

Capacitor  Mass 

Me 

(kg) 

17231.29 

15324.15 

Shield  Mass 

M^jj 

(kg) 

8272.426 

8272.031 

Magnet  Mass 

Mb 

(kg) 

67.55 

67.55 

Total  Mass 

Mtoi 

(kg) 

3,69x105 

3.59x105 

Total  Thmst 

F 

(N) 

4.73x105 

4.84x105 

Thrust-to-Weight 

FAV 

0.131 

0.137 

Specific  Impulse 

_ ki _ 

(s) 

1583.344 

1607.442 

51 


<5) 


niRUSI-TO  -Wr  IfiMI  RATIO 


o.so  n.ia  7.30  11. so 

PRCPSlLLaNT  npSS  PLCIu  RoT^  CKG-'S) 


Figure  27 

Thrust-to-Weight  Ratio  vs.  Propellant  Mass  Flow  Rate 
(Av=5kin/s) 
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Figure  31 

Thrust-tO'Weight  Ratio  vs.  Propellant  Mass  Flow  Rate 
(Av=40kmys) 


CONCLUSIONS  AND  RECOMMENDATIONS 


An  equivalent  circuit  representation  of  the  dense  plasma  focus  device  was 
developed  for  use  in  a  1-D  transient  code  which  solved  for  the  sheath  dynamics  during  the 
rundown  phase  as  well  as  the  current  history  of  the  device.  A  leakage  current  component 
was  inserted  into  the  equivalent  circuit  to  account  for  sheath  current  losses  during  the 
operation  of  the  plasma  focus.  The  leakage  component  was  modeled  as  a  constant 
resistance  which  shunted  part  of  the  external  current  away  from  the  plasma  sheath.  The 
parameter  values  of  the  sheath  were  modeled  as  time  dependent  variables  to  account  for  the 
changing  behavior  of  the  arc  during  rundown.  The  dynamic  nature  of  the  sheath 
parameters  made  it  necessary  to  remodel  the  circuit  after  updating  the  sheath  inductance  for 
each  time  step.  Relations  for  the  plasma  sheath  parameters  were  developed  as  well  as  a 
snowplow  model  for  predicting  the  plasma  sheath  dynamics  during  rundown.  These 
relations  were  inserted  into  the  transient  code  and  another  test  case  was  run  using  the 
Livermore-I  plasma  focus  experiment  for  comparison. 

The  transient  code  used  to  calculate  the  current  history  for  a  given  electrode 
geometry  utilized  an  equivalent  circuit  representation  of  the  DPF  coupled  with  plasma 
relations  for  the  sheath.  This  scheme  solved  the  system  of  coupled  circuit  equations 
numerically  using  LU  decomposition.  The  transient  code  was  tested  on  a  trial  circuit  with 
static  resistance  and  reactance.  These  results  compared  very  well  to  results  calculated 
analytically  and  with  the  SPICE  circuit  modeling  program.  The  error  between  the  transient 
code  calculated  solutions  and  the  analytical  and  SPICE  solutions  were  generally  on  the 
order  of  0. 1  percent  The  accuracy  of  the  transient  solver  was  validated  for  these  series  of 
static  tests.  The  calculated  results  using  the  transient  code  predicted  an  external  current 
history  that  was  in  good  agreement  with  the  experimental  results.  The  calculated  leakage 
current  was  also  in  good  agreement  with  the  experimental  data  obtained  from  the  Livermore 
I  experimental  report  This  experimental  leakage  component  is  an  nf erred  value  that  could 
not  be  directly  measured  in  the  Livermore  I  experiment.  The  presence  of  the  leakage 
current  is  needed  in  this  equivalent  circuit  model  to  account  for  a  current  loss  mechanism 
which  degrades  the  current  delivered  to  the  plasma  sheath.  This  leakage  component  is 
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assumed  to  be  a  leakage  current  over  the  insulator,  but  further  experimentation  is  needed  to 
substantiate  this  assumption. 

After  the  static  circuit  validation  tests  were  conducted,  the  crowbar  branch  of  the 
circuit  was  switched  in  to  determine  its  effect  on  the  load  current  history.  The  results  of 
this  test  exhibited  a  positive  effect  on  the  reduction  of  current  damping  in  the  load.  The 
degree  of  reduction  of  the  current  damping  effect  depends  heavily  on  the  choice  of  crowbar 
inductance  and  resistance  values.  Theoretically,  it  is  possible  to  greatly  reduce  the  damping 
effect  after  peak  current  has  been  reached  by  reducing  the  crowbar  inductance.  However, 
there  exists  a  practical  limit  as  to  how  small  an  inductance  can  realistically  be  achievable  for 
the  crowbar  branch. 

The  Livermore-I  geometry  was  taken  and  modified  in  order  to  provide  the  very  high 
current  that  was  needed  to  heat  the  pinch.  Various  radial  and  axial  anode  variations  were 
tested  in  order  to  determine  the  effect  on  current  history.  It  was  discovered  that  the 
Livermore-I  experimental  configuration  would  not  provide  the  necessary  current  required  to 
produce  fusion  ignition.  As  a  result,  the  charging  voltage  had  to  be  increased  in  order  to 
provide  a  greater  current  delivery  capability.  The  dimensions  of  the  elecU'odes  were  varied 
for  each  increased  value  of  charging  voltage.  These  data  were  recorded  and  the  optimal 
configurations  were  chosen  for  future  testing.  These  configurations  were  chosen  for 
minimal  charging  voltage  requirements  and  minimized  electrode  dimensions.  One 
enhanced  electrode  and  charging  voltage  combination  was  chosen  for  each  operating 
current  requirement  (10,  15,  and  20MA).  Each  of  these  combinations  were  run  on  the 
transient  code  and  the  results  were  input  into  a  DPF  system  code  which  calculated  the 
performance  of  the  thruster.  For  the  best  case  scenario  (lowest  charging  voltage  324  kV, 
highest  sheath  current  20  MA)  the  capacitor  energy  discharge  curve  was  plotted  (Figure 
A.  17).  For  this  case,  it  can  be  seen  that  approximately  2/3  of  the  initial  energy  contained  in 
the  capacitor  is  discharged  in  the  rundown  phase.  A  fraction  of  this  quantity  is  deposited 
into  the  plasma  when  the  pinch  is  formed,  but  a  more  detailed  modeling  of  the  pinch  region 
is  needed  to  determine  this  quantity. 

The  goal  of  this  study  was  to  determine  the  effects  of  changing  electrode  geometry 
on  current  history  and  plasma  focus  performance  in  the  10-20  MA  range.  It  was  found  that 
the  electrodes  did  not  have  to  undergo  drastic  changes  in  configuration  to  achieve  these 
high  currents,  but  the  charging  voltage  had  to  be  increased  to  a  level  which  is  beyond  what 
is  currently  achievable.  The  SHIVA  experiment  conducted  at  Kirtland  Air  Force  Base 
claims  a  load  current  of  approximately  9  MA  at  a  charging  voltage  of  125  kV[10].  This 
provides  a  baseline  as  to  the  current  state  of  high  current  pulse  power  technology.  The 
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DPF  system  would  need  at  least  a  25%  higher  charging  voltage  and  a  discharge  circuit  that 
could  handle  the  high  current  effects.  If  such  a  circuit  could  be  designed  and  implemented, 
perhaps  experimentation  of  the  dense  plasma  focus  in  this  current  range  could  be 
conducted. 

The  effective  result  of  the  new  enhanced  electrode  configurations  is  a  decrease  in 
the  necessary  amount  of  capacitor  mass.  This  allows  the  thrust-to-weighi  ratio  to  be 
increased  significantly  from  previous  calculations.  Another  effect  of  the  enhanced  electrode 
design  is  an  increase  in  fusion  power  provided  by  the  device.  A  longer  annular  region 
allows  more  fill  gas  to  be  entrained  during  rundown.  This  increases  the  amount  of  gas 
entrained  in  the  high  temperature  pinch  region,  which  results  in  a  greater  amount  of  fusion 
fuel  being  burned.  The  thrust-to-weight  ratios  are  maximized  for  short  mission  Av 
requirements.  This  increase  in  thrust-to-weight  ratio  will  decrease  trip  time  of  the  mission 
vehicle  and  its  occupants.  Specific  impulse  for  the  enhanced  electrode  configuration  differs 
little  from  the  Livermore-I  values,  though  for  the  20MA  case,  the  enhanced  configuration  is 
superior  to  the  Livermore-I  configuration. 

Tests  utilizing  different  anode  tip  geometries  were  also  run  using  the  transient 
modeling  code.  The  three  different  anode  tip  geometries  that  were  tested  provided  a  basis 
for  observing  the  effects  of  changing  anode  tip  configurations.  Equilateral,  cylindrical,  and 
extended  triangular  anode  tips  were  tested  for  the  resulting  current  histories  in  the  MA 
range.  Illustrations  of  these  tip  geometries  are  in  Appendix  A  along  with  the  resulting 
current  histories  and  performance  plots.  The  equilateral  and  cylindrical  tips  proved  to  be 
nearly  identical  in  performance  to  the  Livermore-I  sloped  anode  tip,  while  the  extended 
triangular  tip  appeared  to  provide  degraded  performance  characteristics.  The  lack  of  current 
capacity  of  the  extended  triangular  tip  appears  to  be  due  to  the  increased  plasma  impedance 
caused  by  an  elongation  of  the  plasma  sheath  as  it  travels  towards  the  anode  tip.  In  contrast 
with  this,  the  equilateral  and  cylindrical  tips  would  not  experience  the  same  effect  due  to  the 
shorter  lengths  of  the  anode  for  these  geometries.  One  must  note  that  these  conclusions  are 
only  valid  when  analyzing  the  current  carrying  histories  of  the  different  geometries. 
Perhaps  the  extended  tip  geometry  would  be  more  favorable  to  pinch  formation  than  either 
the  equilateral  or  cylindrical  tips. 

The  dense  plasma  focus  device  has  the  possibility  of  becoming  a  desirable  space 
propulsion  concept.  If  the  current  and  capacitor  scaling  laws  hold,  and  if  tougher  higher 
temperature  conductors  can  be  developed,  the  DPF  can  provide  system  performance  that 
will  exceed  other  alternative  concepts.  Questions  about  the  scaling  laws  and  the  physics  of 
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the  plasma  pinch  phase  must  be  studied  further  before  the  DPF  can  be  considered  as  a 
plausible  fusion  energy  source. 

The  results  of  this  study  provided  a  more  realistic  view  of  the  electrode  and  input 
requirements  needed  to  produce  very  high  currents  during  rundown.  These  modifications 
produced  a  more  optimistic  view  of  DPF  thruster  performance  for  different  Av 
requirements.  The  significant  reduction  in  capacitor  mass  is  a  result  of  the  geometrical 
dependence  of  the  capacitor  mass  scaling  law.  This  mass  reduction  plays  a  key  role  in  the 
boosting  of  the  system  thrust-to-weight  ratio.  Further  experimentation  is  needed  to 
substantiate  the  assumptions  that  the  sheath  current  will  not  saturate  in  the  10-20  MA  range. 
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APPENDIX  A 


Various  Electrode  Configurations  and  Associated  Performance  Parameters. 
Figures  A2  through  A.4  which  were  illustrated  earlier  in  the  main  text  are 
are  shown  here  in  repetition.  This  is  done  in  order  to  compare  them  with  the 
different  anode  geometries  presented  in  this  Appendix. 
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Figure  A.  1 

Livermore-I  Plasma  Focus  Electrode  Geometry 
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Figure  A.3 

Thrust-to-Weight  Ratios  for  Enhanced  and 
Livermore-I  Electrodes  (Av=5km/s)  (Same  as  Figure  27) 
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Cathode 
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Anode  with  Equilateral  Tip 


Figure  A.5 

Equilateral  Anode  Tip 
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Figure  A.7 

Thrust-to-Weight  Ratios  for  Equilateral  Anode  Tip 
and  Livermore-I  Electrode  (Av=51an/s) 
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SPECinC  IMPULSE  (S) 


Cathode 


Anode  with  Cylindrical  Tip 


Figure  A.9 

Cylindrical  Anode  Tip 
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Figure  A.  1 1 

Thrust-tO'Weight  Ratios  for  Cylindrical  Anode  Tip 
and  Livennore-I  Electrode  (Av=51cm/s) 
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Figure  A.  12 

Specific  Impulse  for  Cylindrical  Anode  Tip 
and  Livermore-I  Electrode  (Av=5km/s) 


Anode  with  Extended  Triangular  Tip 


Figure  A.  13 

Extended  Triangular  Anode  Tip 
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Figure  A.  15 

Thrust-to-Weight  Ratios  for  Extended  Triangular  Anode 
and  Livermore-I  Electrode  (Av=5km/s) 
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APPENDIX  B 


A  listing  of  the  program  used  to  calculate  the  current  history  and  sheath  parameters 
for  the  rundown  phase  of  plasma  focus  operation  is  included  in  the  following  pages.  The 
code  TRAN.f  is  written  in  FORTRAN  77  and  utilizes  a  transient  circuit  solver  to  calculate 
the  voltage  and  current  histories  for  the  equivalent  circuit  The  plasma  circuit  parameters 
are  calculated  using  suitable  physical  models  in  the  subroutine  SOLVE.  While  the 
rundown  velocity  is  calculated  using  the  snowplow  model. 

Operating  Instructions  for  Using  the  Transient  Code  TRAN.f 

1.  Compile  TRAN.f  using  command  "f77  TRAN.f 

2.  Change  input  parameters  in  input  deck  TRAN.in 

3.  Execute  program  by  simply  typing  "a.ouf ,  program  will  automatically  read  input  data 

and  output  data  into  files: 

i)  lEX.OUT  -  This  file  contains  the  external  circuit  current 

ii)  EPLAS.OUT  -  This  fde  contains  the  plasma  sheath  current 

iii)  IND.OUT  -  This  file  contains  the  plasma  sheath  inductance  history. 

iv)  VOLT  OUT  -  This  file  contains  the  node  voltage  over  the  plasma  sheath. 

v)  VRUN.OUT  -  This  file  contains  the  rundown  velocity  for  the  plasma  sheath. 
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c*  Program  TRAN.F 

c*  Programmer;  Glen  T.  Nakafuji 

c*  Purpose: 

c*  1.  Calculate  current  and  voltage  responses  for  equivalent  circuit  representation  of 
c*  the  dense  plasma  focus  device. 

c*  2.  Calculate  dynamic  circuit  parameter  values  for  plasma  sheath  during  the  rundown 
c*  phase  of  operation. 

c*  3.  Calculates  rundown  velocity  of  propagating  plasma  sheath 
c*  Models: 

c*  1.  Companion  circuit  model  used  to  represent  reactive  elements  in  equivalent  circuit 
c*  2.  Snowplow  model  used  to  calculate  rundown  velocity  of  a  totally  absorbing  arc 
c*  sheet  as  it  propagates  down  the  annular  region. 

c*  3.  Dynamic  models  for  sheath  inductance,  sheath  resistance 


program  TRAN 

dimension  RI{20),Iin(20),LI(20),Vin(2()) 
dimension  Y(20,20),I(20),V(20).Rad(4).ZA{2) 
double  precision  RI,Iin,LI,Vin,Y,I,2A 
double  precision  V,h,c,time,2,Rad,dz,tau 
double  precision  pi,rhoi,mu,imass 

PARAMETER  {pi=3. 1 4 1 59265,mu=l  .256637E-6, vmax=3.5E5) 
COMMON  /  /rhoi.imass 
integer  flag,D,iier 

c*****  Closes  old  output  file 

open  (unit  =  30,file  =  ’scurr.out', status  =  'unknown') 
close  (30,status  =  'delete') 

open  (unit  =  40,file  =  'Icurr.out'.status  =  'unknown') 
close  (40,status  =  'delete') 

open  (unit  =  50, file  =  'volt.out’.status  =  'unknown') 
close  (50,status  =  'delete') 

open  (unit  =  70, file  =  'vrun. out' .status  =  'unknown') 
close  (70,status  =  'delete') 

open  (unit  =  90,file  =  'ind.out'.statys  =  'unknown') 
close  (90,status  =  ’delete’) 

open  (unit  =  100,file  =  'N.oui',status  =  'unknown') 
close  (UX),status  =  'delete') 

open  (unit  =  200,file  =  'convert.out'.status  =  'unknown') 
close  (20G,status  =  'delete') 
c*****  Open  TEMPORARY  input  file 

open  (unit  =  20 .file  =  ’v4.in',status  =  'old') 
c*****  Open  TEMPORARY  output  file 

open  (unit  =  30,file  =  'scurr.out' .status  =  'new') 
open  (unit  =  40,file  =  'Icurr.out'.status  =  'new'i 
open  (unit  =  SO.file  =  'voll.out',status  =  'new') 
open  (unit  =  70,file  =  'vrun.out' .status  =  'new') 
open  (unit  =  90,file  =  'ind.out' .status  =  'new') 
open  (unit  =  lOO.file  =  'N.out' .status  =  'new') 
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open  (unit  =  2CK),file  =  'convert.oui'.status  =  'new') 

c*****  Set  tlag=3  for  Imax  check  prior  to  main  problem 
flag=:3 

c*****  Set  iteration  counter  to  zero 
iter  =  0 

c*****  Set  time  constant  tau 
tau  =  0.0 

c*****  Read  input  data  from  deck 

call  INPUT(RI,Iin,LI.Vin,h,c.Rad,dz,ZA) 


c*****  Initialize  time  and  displacement 
call  INITl(time.z,I,Y,V,dz) 
write  (6,900)  time,z,flag 
10  flag  =  3 

c*****  Initialize  conductance  matrices  Y1  and  Y2  and  I 
call  INIT2(RI,h,c,Vin,Iin,LI,Y.I.flag.D) 
c*****  Solve  I  =  YV  system 

caU  LU(D,Y,I,V,time,h) 
c*******Write  data 

write(30,1000)  time,V(7) 
write(40,1000)  time.(V(l)*RI(3)) 
write(50,1000)  time,V(l) 
write(90,l(X)0)  time,LI(3) 
c******  INPUT  check  mode  only 
c  call  CHECK(V,flag,lime,tau) 

call  SOLVE(I,V,time,RI,U,h,z,ZA.Rad,flag,Iin,Vin.dz) 
c*****  Update  screen  every  n  iterations 
c  write  (6,1  (X)0)  time, z 

c****"*  Advance  iteration  counter  by  1 
iter  =  iter  +  I 

c*****  Check  if  endtime  reached 
if  (z.GE.ZA(l))  then 
call  VOL(ZA,Rad,V,iter) 
write(6,*)  iter 
stop 
endif 

c*****  Clear  Arrays  before  next  p;iss 
call  ZERO(I,V,Y,D) 
goto  10 

800  formatCdimension  =',I2) 

900  formai('time  =',E13.7E2,'position  ='.E13.7E2,3x,I4) 
1000  format(E13.7E2,3x,E13.7E2) 
stop 
end 


c  subroutine  INPUT 
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c  Purpose:  Initializes  parameter  arrays  with  initial  values 

subroutine  INPUT(RI,Iin,LI,Vin,h,c,Rad,dz,ZA) 

dimension  RI(20),Iin{20),LI(20),Rad(4) 

dimension  Vin(20XZA(2) 

double  precision  h,c,ZA,RI,Iin,LI,Rad,dz 

double  precision  rl,r2,r3,Vin,pi,mu,vmax,rhoi,imass,area 

PARAMETER  (pi=3.14159265,mu=1.256637E-6.vmax=3.5E5) 

COMMON  /  /rhoi.imass 

real  div 

c****  reads  parameters  from  input  deck 
c****  read  in  1/ro  (initial  conductance) 
read  (20,*) 
read  (20,*)  rl 
RI(1)  =  l./rl 

C++**  jjj  (initial  crowbar  conductance) 
read  (20,*) 
read  (20,*)  r2 
RI(2)  =  l./r2 

C++++  jcad  in  i/rl  (leakage  conductance) 
read  (20,*) 
read  (20,*)  r3 
RI(3)  =  L/r3 

C++++  yj  lq  (initial  circuit  inductance) 
read  (20,*) 
read  (20,*)  LI(1) 

C++++  read  in  Ix  (crowbar  inductance) 
read  (20,*) 
read  (20,*)  LI(2) 

C++++  ij^  initial  capacitor  voltage 
read  (20,*) 
read  (20,*)  Vin(3) 

C++++  ij^  cathode  radius 
read  (20,*) 
read  (20,*)  Rad(l) 

C++++  read  in  anode  radius 
read  (20,*) 
read  (20,*)  Rad(2) 

C++++  j.g^j  jjj  jjQgg  radius 
read  (20,*) 
read  (20,*)  Rad(3) 

C++++  jjj  jijj^g  5|gp 

read  (20,*) 
read  (20,*)  h 

C++++  capacitance  value 

read  (20,*) 
read  (20,*)  c 

C++++  if,  anode  length 
read  (20,*) 
read  (20,*)ZA(1) 

C++**  pgmi  if,  curve  in  length 
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read  (20  ♦) 
read  (20,*)  ZA(2) 
c****  read  in  sheath  thickness 
read  (20,*) 
read  (20,*)  dz 

c****  read  in  coefficient  for  limiting  leakage  current 
read  (20,*) 
read  (20,*)  RI(20) 
c***  read  in  initial  fill  gas  density 
read(20,*) 
read(20,*)  rhoi 

c****  calculate  initial  sheath  inductance 
div  =  Rad(l)/Rad(2) 

LI(3)  =  (2.E-7)*ALOG(div)*(dz  +  I. *.14) 
c****  set  for  now  initial  plasma  resistance 
RI(4)  =  l./l.E-i4 

area=  pi*((Rad(l)**2)  -  (Rad(2)**2)) 
imass  =  rhoi*dz*aiea 

return 

end 

^  ^  :ii  3)1 ^c,|c  ^  *>)<<*<  >K  >0  I#  >*1  *  Ik  *******  t  ** ’*■***  * 

c*  Subroutine  INITl 

c*  Purpose:  initialize  time,  displacement  and  arrays 

^****  *********  *********************;|i;(|,),;t„),,|,,),,),,(,,|,,(,,^;^,t„),:(,,|C,^,(,,),,),,(,;^,|,,^,^,^,^,^^,,^,^ 

subroutine  INITl  (time,z,I,Y,V,dz) 
dimension  I(20),Y(20,20),V(20) 
double  precision  I,Y,V,time,z,dz 
integer  j,k 

c*****  Initialize  displacement  to  1st  thickness  of  sheath 
z  =  dz 

c*****  Zero  out  arrays  and  variables 
time  =  0.0 
do  10j=  1,20 

ia)=o.o 
Va)=0.0 
do  5  k  =  1,20 
Ya.k)  =  0.0 
5  continue 
10  continue 

return 

end 

j,*****************  *********  + 

c*  Subroutine  INIT2 

c*  Purpose:  Initialize  conductance  matrices  and  select  which 

c*  circuit  to  implement  based  on  the  flags  that  are  set. 

c+ 

^**********************************:|l,|t,tl,)ia)t***************i|c^l:|,;^:(,:^,|c****** 
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subroutine  INIT2(RI,h,c,Vin,Iin,LI, Y ,I,flag,D) 
dimension  RI(20),Vin(20),Iin(2{)),LI(20) 
dimension  Y(20,20)J(20),Y1  (20.20), Y2(20,20) 

double  precision  RI,Vin,Iin,LI,Y,I,h,c 
integer  flag,D,j,k 

if  (flag.NE.2)  then 

c******  Initialize  elements  of  7x7  matrix  for  circuit  w/o  crowbar 
Y1(1,1)  =  RI(1)  +  RI(3) 

Yl(l,2)=-RI(l) 

Yl(l,7)=  1. 

Yl(2,l)  =  -RI(l) 

Yl(2,2)  =  RI(1) 

Yl(2,6)  =  -1. 

Yl(3,5)=  1. 

Yl(3,6)  =  1. 

Y1(4,4)  =  RI(4) 

Yl(4.7)  =  -l. 

Yl(5,3)  =  (2*c)/h 
Yl(5,5)  =  -1. 

Yl(6.2)  =  -h/(2*LI(l)) 

Yl(6,3)  =  h/(2*LI(1)) 

Yl(6,6)  =  -1. 

Yl(7,l)  =  h/(2*LI(3)) 

Yl(7,4)  = -h/(2*LI(3)) 

Yl(7,7)  =  .l. 

c*****  Initialize  and  update  solution  vector  I 
1(5)  =  lin(l)  +  ((2*c)/h)*Vin(3) 
cl=-h/(2*LI(l)) 

1(6)  =  cl*(Vin(3)  -  Vin(2))  -  Iin(2) 
c2  =  -h/(2*LI(3)) 

1(7)  =c2*(Vin(l)  -  Vin(4))  -  Iin(4) 

D  =  7 

c*****  Write  Y1  into  general  Y  conductance  matrix 
do  200j  =  l.D 
do200k=  l,D 
Y(j.k)  =  Yl(jJt) 

200  continue 

endif 

if  (flag.EQ.2)  then 

c*****  Initialize  Y2  matrix  if  flag  =  2 
Y2(1,1)  =  RI(1)  +  RI(3) 

Y2(1,2)  =  -RI(1) 

Y2(l,8)=  1. 

Y2(l,9)=  1. 

Y2(2,1)  =  -RI(1) 

Y2(2.2)  =  RI(1) 
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Y2(2,7)  =  -1. 

Y2(7,6)  =  1. 

Y2(7,7)  =  1. 

Y2(4,4)  =  RI{4) 

Y2(4,9)  =  -1. 

Y2(5,5)  =  RI(2) 

Y2(5,8)  =  -1. 

Y2(6,3)  =  (2*c)/h 
Y2(6,6)  =  -l. 

Y2(3,2)  =  -h/(2*LI(l)) 

Y2(3,3)  =  h/(2*LI(l)) 

Y2(3,7)  =  -1. 

Y2(8,l)  =  h/(2*LI(2)) 

Y2(8,5)  =  -h/(2+LI(2)) 

Y2(8,8)  =  -l. 

Y2(9,l)  =  h/(2*LI(3)) 

Y2(9,4)  =  -h/(2*LI(3)) 

Y2(9,9)  =  -1. 

c*****  Initialize  I  vector  tor  Y2  system 
1(6)  =  lin(l)  +  ((2*c)/h)*Vin(3) 
cl  =-h/(2*LI(l)) 

1(3)  =  cl*(Vin(3)  -  Vin(2))  -  Iin(2) 
c3  =  -hy(2*LI(2)) 

1(8)  =  c3*(Vin(l)  -  Vin(5))  -  Iin(3) 
c2  =  -h/(2*LI(3)) 

1(9)  =  c2*(Vin(l)  -  Vin(4))  -  Im(4) 

D  =  9 

c*****  Write  Y2  ii.to  general  Y  conductance  matrix 
do300j  =  1,D 
do  300  k=  1,D 
Y(j,k)  =  Y2(j,k) 

300  continue 


endif 


return 

end 

c*  Subroutine  LU 

c*  Purpose;  Solves  I=YV  system  using  LU  decomposition  and  also 
c*  advances  Lime  step 


subroutine  LU(row,A,B,X,time,h) 
dimension  L(20,20),U(20,20),B(20),A(20,20) 
dimension  BP(20),X(20) 
double  precision  L,U,B,A,BP,X,time.h 
integer  row, col 
col=row+l 
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return 

end 

c*  Subroutine  forward 

c*  Purpose;  Part  of  subroutine  LU  to  solve  system 

subroutine  forward(row,L.BP.B) 
dimension  L(20,20),BP(2()).B(20) 
double  precision  L,BP,B 
integer  i,j,row 

BP(1}=B(1)/L(l.l) 
do  200  i=2.row 
temp=0.0 
do  I00j=l,i-1 
temp=temp+L(i,j)*BP(j) 

1(X)  continue 

BP(i)=(B(i)-temp)/L(i,i) 

200  condnue 

return 

end 

c*  Subroutine  back 

c*  Purpose:  Part  of  subroudne  LU  to  solve  system 

subroudne  back(row,U,X,BP) 
dimension  U(20.20).X(20),BP(20) 
double  precision  U,X.BP,sum 
integer  ij.row 

X(row)=BP(row)AJ(row,row) 
do  200  i  =  row- 1,1,- 1 
sum=0.0 

do  180  j=i-f-l,row 
sum  =  sum-t-U(i,j)*X(j) 

180  continue 

X(i)=(BP(i)-sum)/U(i,i) 

200  condnue 


return 

end 


Q**********  ♦****♦  +  +  *  +  *  +  *  ♦**  +  *****  +  **  +  :tt  +  +  +  ****  +  ^*  +  *,(,  +  *****;(,  +  **,),***,(,**,(.,(,^ 

c*  Subroutine  pivot 

Purpose:  Part  of  subroudne  LU  to  solve  system 


subroutine  pivot(row,A,B) 
dimension  A(20,20),B(20) 
double  precision  A,B,check,sub,subl,temp 
integer  row 
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I  temp=0.0 
sub=0.0 
sub  I  =0.0 

check  =  A(row,row) 
if  (check.EQ.O)  then 
do  5  i  =  1  .row 
sub  =  A(l,i) 

A(l,i)  =  A(row,i) 

A(row,i)  =  sub 

5  continue 

subl  =B(1) 

B(l)  =  B(row) 

B(row)  =  subl 
endif 

do  20  i  =  l.row-1 
check  =  A(i,i) 

if  (check,EQ.O)  then 
do  10  j  =  l.row 
sub  =  A(i+ 1  ,j) 

A(i+l,j)  =  A(i.j) 

A(i,p  =  sub 

10  continue 

subl  =  B(i) 

B(i)  =  B(i+l) 

B(i+I)  =  subl 
endif 

20  continue 

c*****  check  diagonal  elements  for  zeros 
do  30  i  =  I  .row 
temp  =  A(i.i) 

c*****  if  zero  go  back  to  beginning  and  swap  rows 
if  (temp.EQ.O)  then 
goto  1 
endif 

30  continue 


return 

end 

:if(  ifc  3)c  :4c  4[  ife  ifc  34c  jfc  ^  ^  ift  4c  3^  *  ^  ^  4t  %  )(i  ^  ^  ^  ^  iti  3(1 

c*  Subroutine  ZERO 

c*  Purpose:  Zero  out  general  matrix  and  vectors  before  next 
c*  iteration 

Q’^**^*^***^^*****^^*^‘^‘**^’¥***^^**************  ************’¥*  ******<¥**** 

subroutine  ZERO(I,V,Y,D) 
dimension  V(20),I(20),Y(20,20) 
double  precision  V.I.Y 
integer  D 
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do  2t)j=l,D 
do  10  k=l.D 
Y(j.k)  =  0.0 
10  continue 

l(j)=0.0 
V(i)=0.0 
20  continue 

return 

end 

c*  Subroutine  CHECK 

c*  Purpose;  Checks  if  current  of  initial  circuit  (7x7)  is 
c^  at  Imax,  if  it  is,  switch  in  9x9  by  setting  flag  to  2 
c*  otherwise  keep  flag  at  1  or  3 

^^t^ti^ifi^iiti*****************************^****************  ************ 

subroutine  CHECK(V, flag. time, tau) 
dimension  V(20) 

double  precision  V.time.tau.temp.max.unax 
integer  flag 

c******  If  flag=3  program  in  Imax  check  mode 
if  (flag.EQ.B)  then 

c******  If  current  value  of  Ipf  is  greater  than  previous  max.  store 
temp  =  ABS(V(7)) 
if  (temp.GT.max)  then 
max  =  temp 
tmax  =  time 
endif 

c*****  Otherwise,  check  for  end  of  discharge,  if  at  end  set  flag  to  1 

if  (time.GE.tau)  then 
flag  =  1 
time  =  0.0 

write  (6,1000)  max,tmax 
endif 

endif 

c*****If  flag=I  program  starts  problem  with  7x7  system 
if  (flag.EQ.l)  then 

c*****Check  if  time  >  unax,  if  it  is  then  switch  circuit  by  flag  =  2 
if  (time.GT.tmax)  then 
flag  =  2 
endif 

endif 

1000  formatCmax  current  or,fl3.7,‘occurs  at  t  =  ',fl5.8) 
return 
end 

^^ltii^^^l^7^:^li^^i^lt^i^^fL*^^^L^^*7^***************,^***********************i|t^L^l^^ll^*:^^* 

c*  Subroutine  SOLVE 

c*  Purpose:  Solve  for  sheath  parameters  and  calculate  current 
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c*  displacement,  and  replace  new  values  into  inii  vectors 

subroutine  SOLVE(I,V,time,RI,LI,h.z,ZA.Rad.flag,Iin.Vin,d/,) 

dimension  I(20).V(20).RI(20)XI(20).Rad{4) 

dimension  Iin(20),Vin(20),ZA(2) 

double  precision  I,V,RI,LI,Rad.iime.h.z 

double  precision  vrun,Iin,Vin,denom,pi,mu.vmax 

double  precision  imass,rhoi.area,Inew 

double  precision  Rlo,dz 

double  precision  Ilift,ZA,r,vcorr,delta 

double  precision  intold.  inisum,  intnew 

PARAMETER  (pi=3. 1 4 1 59265,mu=  1 .256637E-6.vmax=3,5E5) 

COMMON  /  /rhoi,imass 

integer  flag 

c****  Set  initial  leakage  conductance 
Rio  =  RI(3) 

c****  Replace  old  currents  and  voltages  with  new  values 
do  10j=1.4 
Vin(j)  =  V(j) 

10  continue 

Iin(l)  =  V(5) 

Iin(2)  =  V(6) 

Iin(4)  =  V(7) 

Inew  =  V(7) 

c******''‘*Don't  need  9x9 
c  if  (flag.EQ.2)  then 

c  do20j=l.5 

c  Vin(j)  =  V(j) 

c20  continue 

c  Iin(l)  =  V(6) 

c  Iin(2)  =  V(7) 

c  Iin(3)  =  V(8) 

c  Iin(4)  =  V(9) 

c  Inew  =  V(9) 

c  endif 


c*****‘'‘*Calculate  rundown  velocity 

c*****  with  Liftoff  current  calculation 
vrun  =  0.0 

c****  Adjusted  current  for  liftoff 

Ilift  =  (Inew**2)  -  (4144.82852**2) 

c*****  CaU  radius  and  velocity  correction  routine 
call  CORRECT(Rad,ZA,r,vcorr,z) 

^  .<1  >|>  ■*<  %  itiiti  41 « >x  1)1  >)■  iti  ****>)•>*  >l>  *  4»li  >)•«  f  4>  ***>•>  I«>  *  >4<  *  4>  4  4"t< 
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c******  Snowplow  model  using  trapezoidal  approximation 
c****'^***CaJculate  LI**2  integral  using  Trapezoidal  approx, 
if  (Ilift.GT.0.0)  then 
intnew  =  (LI(3)*(Ilift))/z 
delta  =  .5*h*(intold  +  intnew) 
intsum  =  intsum  +  delta 
intold  =  intnew 
area  =  pi‘^(Rad(l)**2  -  r**2) 
denom  =  2*(imass+(rhoi*area*z)) 
vrun  =  (vcorr/denom)*intsum 
endif 

Qj(>****ltl+l)lltt**+**  +  *+  +  l)ll(l+*l(i*++i)t*  +  l(i+*+l^*++*^*++*+*+*+**************l*#* 

c******  Limit  sheath  velocity  to  implosion  velocity 
if  (vrun.GT.vmax)  then 
vrun  =  vmax 
endif 

c***’''***Calculate  new  position 
z  =  z  +  vrun*h 

c***  Write  rundown  velocity  to  output 
write  (70,  KKX))  time. vrun 
c******  Calculate  new  sheath  inductance 
div  =  Rad(l)/r 

LI(3)  =  (2E-7)*ALOG(div)*(z+  l.*.14) 

1000  forraat(E13.7E2,4x,E13.7E2) 
return 
end 

^ Iti ^ «  4c  1*1  ^ ^  1|1  :)c ;ti * Itcift Itcifi  1)1 4i # :ti « :ti « >»*>*>« ****««* 4>  * 4>  **  4< >t<4>  I*!  <«•  « 

c*  Subroutine  CORRECT 

c*  Purpose:  Correct  velocity  and  radius  for  sloped  anode  section 
c* 

subroutine  CORRECT(Rad,ZA,r,vcorr,z) 
dimension  Rad(4),ZA(2) 

double  precision  Rad,21A,r,vcorr,z,rboss,ztip,zcheck 
double  precision  rl,r2,zrad,drdz,zcv,iop.bottom,ra 
ztip  =  ZA(l) 
zcv  =  ZA(2) 
rboss  =  Rad(3) 
ra  =  Rad(2) 
zcheck  =  ztip  -  zcv 
if  (zcheck.LT.0.0)  then 
write(6,*)  ('Electrode  Geometry  Invalid') 
stop 
endif 

zrad  =  zcv  +  ra  -  rboss 

c****  If  anode  hasn't  started  curving  in  keep  factors  constant 
vcorr  =  1.0 
r  =  ra 

c****  If  anode  has  started  curving  in,  calculate  velocity  and 
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c****  radius  adjustment 
if  (Z.GE.ZCV)  then 

if  (z.GE.zrad)  then 
r  =  rboss 
vconr  =  1.0 
goto  10 
endif 

top  =  -z  +  zcv 

bottom  =  SQRT(((ra  -  rboss)^*2)  -  ((z  -  zcv)**2)) 
drdz  =  top/bottom 

c****  correction  coefficient  for  axial  velocity  component 
vcorr=  l./SQRT(l+drdz**2) 
rl  =  (ra  -  rboss)**2 
r2  =  (z  -  zcv)* *2 

c***  correction  for  radial  length  component 
r  =  rboss  +  SQRT(rl  -  r2) 
endif 

10  continue 

return 
end 

^  itDft  ^  ]fi  ^  ^  ^  ^  ^  ^  iff  Ifi  ^  ifi  4t  ^  ^  #  4 1#  ^  4t  1#  4c  ^  ^ 
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c  subroutine  VOL 

c*  Purpose;  Calculates  volume  of  annular  region  in  modified  geometry  and  writes 
c*  pertinent  endtime  values  to  output  file  "convertout" 

Ql^;^l|li^:^tl^^^t|l^^c^l^^,lt‘**’^‘****‘l‘******'H‘***i‘***’^‘*’^‘*^l*0**^l‘***’l‘^‘*****'^*  ******  ****** 

subroutine  VOL(ZA,Rad,V,iter) 
dimension  ZA(2),Rad(4),V(20) 

double  precision  ZA,Rad,areal,area2,volume,zinspi,totlen 
double  precision  V,N 
integer  iter 

PARAMETER  (pi=3.14i59265.mu=1.256637E-6,vmax=3.5E5) 
areal  =  pi*(Rad(l)**2  -  Rad(2)**2) 
area2  =  pi*(Rad(l)**2  -  Rad(3)**2) 
c*****  Insulator  length 
zins  =  .14 

c*****  Total  length  of  anode  for  Livermore  I  assumption 
totlen  =  .382 

c*****  Compute  Volume  of  Annulus 
volume  =  areal  *21A(2) 
write  (200,*)  ('Volume  is:’) 
write  (200,10()0)  volume 
c*****  Write  final  current  of  sheath 

write  (200,*)  ('Final  current  in  sheath  is;') 
write  (200,1000)  V(7) 
c*****  Write  number  of  iterations 

write  (200,*)  (Total  number  of  iterations  is:’) 
write  (2(X),20{)0)  iter 
1000  format  (El 3.8E2) 

2000  format  (15) 
return 
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